# Bayesian Physics-Informed Neural Networks for Robust System Identification
of Power Systems

Simon Stock1, Jochen Stiasny2, Davood Babazadeh1, Christian Becker1, Spyros
Chatzivasileiadis2 1Institute of Power and Energy Technology, Hamburg
University of Technology, Hamburg, Germany
2Department of Wind and Energy Systems, Technical University of Denmark, Kgs.
Lyngby, Denmark
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This paper introduces for the first time, to the best of our knowledge, the
Bayesian Physics-Informed Neural Networks for applications in power systems.
Bayesian Physics-Informed Neural Networks (BPINNs) combine the advantages of
Physics-Informed Neural Networks (PINNs), being robust to noise and missing
data, with Bayesian modeling, delivering a confidence measure for their
output. Such a confidence measure can be very valuable for the operation of
safety critical systems, such as power systems, as it offers a degree of
‘trustworthiness’ for the neural network output. This paper applies the BPINNs
for robust identification of the system inertia and damping, using a single
machine infinite bus system as the guiding example. The goal of this paper is
to introduce the concept and explore the strengths and weaknesses of BPINNs
compared to existing methods. We compare BPINNs with the PINNs and the
recently popular method for system identification, SINDy. We find that BPINNs
and PINNs are robust against all noise levels, delivering estimates of the
system inertia and damping with significantly lower error compared to SINDy,
especially as the noise levels increases.

###### Index Terms:

Bayesian Physics-Informed Neural Networks, System Identification, Swing
Equation, Power System Dynamics

##  I Introduction

Power systems dynamic behavior is drastically changing due to the rapid
integration of converter-connected energy sources, and the gradual
decommissioning of conventional power plants. This leads to overall lower
levels of inertia and a less damped system. Additionally, the uncertain and
fluctuating output of renewable energy sources, such as wind farms and solar
PVs, results in frequent changes of the power system dynamic characteristics.
This requires advanced approaches to assess the power system state rapidly,
continuously, and with confidence, in order to determine the current level of
inertia and damping and commit accordingly the necessary frequency reserves to
avoid instability and blackouts. The proliferation of different types of
sensors, such as phasor-measurement units (PMUs), smart meters, and others,
has enabled the high-frequency monitoring of power system behavior through the
collection of significant amounts of data. Such data can help towards the so-
called online system identification, for which multiple approaches exist in
the literature as we attempt to summarize below. To determine the system
dynamics, the temporal derivative of the data is usually calculated. This
helps estimate the parameters of a set of differential equations which
describe the system behavior.

System identification has so far used classical linear approaches such as
Dynamic Mode Decomposition (DMD) [1] or Koopman theory [2], and approaches
that rely on filtering methods such as Kalman filtering [3]; recently, there
has been increased focus on parsimonious approaches, such as SINDy [4] to cope
with nonlinearities. We briefly detail SINDy in Section II-C of this paper.
Nevertheless, SINDy is prone to noise in the data due to its point-wise
differentiation that could lead to large errors. In that respect, Physics-
Informed Neural Networks have appeared as a very promising alternative [5], as
they fit the whole trajectory before differentiation; this helps with
smoothing the values before getting differentiated, which makes them robust
against noise. For example, [5] shows that under strong nonlinearities and
noise, Physics-Informed Neural Networks (PINN) can estimate the inertia and
damping of the system with small error.

However, SINDy and PINNs are by design not able to give a statement about the
confidence level of their estimates. In contrast, Bayesian approaches have the
capability to quantify uncertainty as they try to fit a whole family of
functions instead of a single function. By inferring a posterior distribution
of the parameters from the measurement data, the mean value can be seen as the
best guess of the estimate and the standard deviation as a measure of its
uncertainty. Following this path, [6] has presented an approach for system
identification of power system dynamics utilizing Bayesian inference. However,
this approach requires prior knowledge about the region of the target
parameter.

To the best of our knowledge, this paper introduces for the first time the
Bayesian Physics-Informed Neural Networks for use in power system
applications. Bayesian Physics-Informed Neural Networks (BPINNs) combine the
advantages of the PINNs, being robust to noise and missing data, with Bayesian
modelling that is able to deliver a confidence measure of the estimate,
through the standard deviation of the posterior. In contrast with existing
approaches, BPINN does not require prior knowledge about either the target
values (parameters) and their potential range or the noise in the measurement
data, i.e. fidelity of the sensors [7].

BPINNs can have a wide range of applications in power systems, going well
beyond robust system identification that we explore in this paper. The goal of
this paper is to introduce the concept and explore the strengths and
weaknesses of BPINNs compared to standard PINNs and the powerful conventional
approaches. Enriching, however, Physics-Informed Neural Networks, with a
confidence measure about their estimate – as Bayesian PINNs do – is expected
to open a wide range of opportunities for application in power systems.

This paper is structured as follows. Section II formulates the problem,
provides the related equations, and presents the proposed approach in detail.
Section III investigates the BPINN and compares its performance to SINDy and
the PINN. Section IV concludes.

##  II Methodology

In this section, we introduce the Bayesian Physics-Informed Neural Network
(BPINN) formulation used throughout this paper. To do that, first, we briefly
revisit the neural network and PINN formulations, and then we present the
Bayesian PINN and highlight similarities between the deterministic and
Bayesian formulations.

Assume a dynamic model described by the following differential equations (in
Section II-B, we detail the formulation for the specific power system dynamics
problem):

$\displaystyle\dot{\boldsymbol{x}}$
$\displaystyle=f(\boldsymbol{x},\boldsymbol{u};\boldsymbol{\lambda})$ (1)

with solution $\boldsymbol{x}(t,\boldsymbol{u})$, $\boldsymbol{x}$
representing the states of the system and $\boldsymbol{u}$ the inputs of the
system. $\boldsymbol{\lambda}$ describes the system parameters that are to be
determined and $f$ is the operator mapping the system parameters to the
states. Neural networks are generally used as function approximators for
various problems. In the case of dynamic systems, their goal is to determine a
surrogate model ${g}(t)$ that maps a time-dependent input vector to the target
trajectory of the states $\boldsymbol{x}$ as follows:

$\displaystyle
g(t;\boldsymbol{\Theta})=\hat{\boldsymbol{x}}(t;\boldsymbol{\Theta})\approx\boldsymbol{x}(t,\boldsymbol{u};\boldsymbol{x_{0}},\boldsymbol{\lambda})$
(2)

Here, $\boldsymbol{x_{0}}$ describes the initial state of the system. The
surrogate model should be determined using a set of measurement data
$\mathcal{D}$ with
$\mathcal{D}=\\{(\boldsymbol{x}^{(i)},\boldsymbol{u}^{(i)})\\}_{i=1}^{N}$.
During training, the distance between $\hat{\boldsymbol{x}}(t)$ and the target
trajectory is minimized by varying the neural network parameters
$\boldsymbol{\Theta}$, i.e. the NN weights and biases. Equation (3) presents
the formulation of such an optimization problem, using the root-mean-squared
error as the distance measure. The training data
$\boldsymbol{x}\in\mathcal{D}$ follow an additive Gaussian noise model
$\boldsymbol{\eta}$ with unknown standard deviation. Such a Neural Network
training procedure is a supervised learning problem.

$\displaystyle\min_{\Theta}\frac{1}{N}\sum_{i=1}^{N}\sqrt{(\hat{\boldsymbol{x}}^{(i)}(\Theta)-(\boldsymbol{x}_{true}^{(i)}+\boldsymbol{\eta}^{(i)}))^{2}}$
(3)

Physics-Informed Neural Networks (PINNs) have been proposed to determine the
parameters $\boldsymbol{\lambda}$ [8] utilizing automatic differentiation to
find $\frac{d}{dt}\hat{\boldsymbol{x}}$ and feed this into a physics
regularization term:

$\displaystyle h(t,\boldsymbol{u};\boldsymbol{\Theta},\boldsymbol{\lambda})$
$\displaystyle=\frac{d}{dt}\hat{\boldsymbol{x}}-f(\hat{\boldsymbol{x}},\boldsymbol{u};\boldsymbol{\lambda})\stackrel{{\scriptstyle!}}{{=}}0$
(4)

The goal is to determine the parameters $\boldsymbol{\lambda}$ of a known
function $f$. This, in theory, does not require any external data, and,
therefore makes the problem an unsupervised learning problem, as shown in [5].
The Physics-Informed Neural Networks usually combine (3) and (4) in the loss
function of the neural network training, to determine the system states
$\boldsymbol{x}$ and parameters $\boldsymbol{\lambda}$.

###  II-A Bayesian PINN

Bayesian Neural Networks (BNNs) were introduced about three decades ago to
address the problem of uncertainty in data, and were shown to deliver good
estimates in the presence of noise [9]. In BNNs, the weights, biases and
outputs of the neural network are formulated as probability distributions and
during the training phase, Bayes’ theorem is used to infer their posterior
distributions $P(\boldsymbol{\Theta}|\mathcal{D})$ based on the data
$\mathcal{D}$ and the specified prior $P(\boldsymbol{\Theta})$ of the
parameters. As a result, the Bayesian neural network can be seen not as a
single function, but as a family of functions that also deliver a confidence
measure about their function approximation based on the variance of the
posterior. This formulation also leads to a different loss formulation
compared to (3) which is now based on the following equation.

$\displaystyle
P(\mathcal{D}|\boldsymbol{\Theta})=\prod_{i}^{N}\frac{1}{\sqrt{2\pi\sigma_{x}^{(i)^{2}}}}\text{exp}(-\frac{(\hat{\boldsymbol{x}}^{(i)}(\boldsymbol{\Theta})-\boldsymbol{x}^{(i)})^{2}}{2\sigma_{x}^{(i)^{2}}})$
(5)

To keep the approach widely applicable, here we do not assume the measurement
error to be known a priori; this means that $\sigma_{x}$ is also a parameter
that has to be determined during training. The posterior distribution of
parameters $\boldsymbol{\Theta}$ can then be calculated based on

$\displaystyle
P(\boldsymbol{\Theta|}\mathcal{D})=\frac{P(\mathcal{D}\boldsymbol{|\Theta})P(\boldsymbol{\Theta})}{P(\mathcal{D})}$
(6)

using a Bayesian inference algorithm. Still, this formulation is only
applicable to forward problems where the goal is to approximate the temporal
sequence of states $\boldsymbol{x}(t,\boldsymbol{u})$. In a power systems
context, this allows us to estimate e.g. the evolution of frequency or rotor
angle, but it does not allow us to perform a system identification. To
determine the parameters $\boldsymbol{\lambda}$, we need to introduce a
physical regularization similar to (4), as proposed in [7] for Bayesian
approaches, introducing the so-called Bayesian Physics-Informed Neural Network
(BPINN). By doing that, not only has the trajectory estimate from the BNN to
fit the true trajectory, but also has to validate the physical equation. As a
result, the loss function for the BPINN training becomes [notice (7) includes
(5)]:

$\displaystyle
P(\mathcal{D}|\boldsymbol{\Theta},\boldsymbol{\lambda})=P(\mathcal{D}|\boldsymbol{\Theta})\prod_{i}^{N}\frac{1}{\sqrt{2\pi\sigma_{h}^{(i)^{2}}}}\text{exp}(-\frac{(h^{(i)}(\boldsymbol{\Theta},\boldsymbol{\lambda}))^{2}}{2\sigma_{h}^{(i)^{2}}})$
(7)

The joint posterior of $\boldsymbol{\Theta}$ and $\boldsymbol{\lambda}$ can be
determined following Bayes’ theorem based on the prior distributions of
$P(\boldsymbol{\Theta)}$ and $P(\boldsymbol{\lambda)}$.

$\displaystyle
P(\boldsymbol{\Theta,\lambda|}\mathcal{D})=\frac{P(\mathcal{D}\boldsymbol{|\Theta,\lambda})P(\boldsymbol{\Theta,\lambda})}{P(\mathcal{D})}$
(8)

This process of a Bayesian Physics-Informed Neural Network is illustrated in
Fig. 1.

Figure 1: Bayesian Physics-Informed Neural Network schematic

However, the posterior distribution of unknown parameters
$P(\boldsymbol{\Theta,\lambda|}\mathcal{D})$ is intractable. To find it, we
use variational inference, which formulates this problem as an optimization,
so computationally efficient optimization libraries can be used. To solve this
optimization problem, the Stein Variational Gradient Descend (SVGD) algorithm
is used in this paper as proposed in [10].

###  II-B Simulation

To help our understanding about the BPINN performance, we introduce the
Bayesian Physics-Informed Neural Networks for power systems using a single-
machine infinite-bus system, as follows:

$\displaystyle\dot{\delta}$ $\displaystyle=\Delta\omega$ (9)
$\displaystyle\dot{\Delta\omega}$
$\displaystyle=\frac{1}{m}(P-d\,\Delta\omega-B\,\text{sin}\delta)$ (10)

with $m$ being the inertia, $d$ the damping, $B$ the entry of the susceptance
matrix and $P$ the active power. The states of the system are the angle and
the frequency deviation: $\boldsymbol{x}=\\{\delta,\Delta\omega\\}$. A single-
machine infinite bus system can also often represent an aggregate of a larger
power system; estimating the parameters $\boldsymbol{\lambda}=\\{m,d\\}$
through the BPINN delivers an estimate for the system inertia and damping,
including a confidence measure about the estimate. We studied four different
scenarios ranging from fast system dynamics to slow system dynamics as
described in Table I. The corresponding trajectories, assuming no measurement
noise, are shown in Fig. 2. For all scenarios, $B=0.2$ and $P=0.1$, while the
initial states are $\boldsymbol{x}_{0}=\\{0,0\\}$.

TABLE I: Evaluation scenarios Scenario |  $m$ in p.u. |  $d$ in p.u.
---|---|---
fast dynamics 1 (fd1) | $0.3$ | $0.15$
fast dynamics 2 (fd2) | $0.6$ | $0.3$
slow dynamics 1 (sd1) | $1.4$ | $1.1$
slow dynamics 2 (sd2) | $1.7$ | $1.4$
Figure 2: $\omega$ and $\delta$ trajectories of all scenarios ($0\,\%$ noise)

For all tests, a BPINN with 10 neurons and 1 hidden layer, and a standard
trajectory length of $T=27\,s$ was used at a sampling frequency of
$10\,\text{Hz}$. The inputs of the BPINN are $t$ and $P$. To avoid an
initialization bias, 10 runs are performed for all tests and the average mean
and average standard deviation are used for the BPINN results presented in
Table II and Table III. For SINDy and the PINN 10 runs are performed per
scenario and the average of the estimates is presented in the tables Table II,
Table III. In contrast to BPINNs, which also output the standard deviation of
their estimate in each run, notice that SINDy and PINN do not offer any
confidence measure, they only offer a point estimate. The BPINN was
implemented in Python using packages $pytorch$ and $numpyro$; the code can be
found in GitHub [11].

###  II-C SINDy algorithm

We compare our approach with one of the recently popular approaches for system
identification: the Sparse Identification of Nonlinear Dynamics (SINDy)
algorithm, proposed in [4]. This approach is part of the parsimonious
approaches which focus on the active terms of the dynamic systems, to reduce
the computational effort and produce explainable results. To determine a
generalized model with as few active terms as possible, sparse regression is
utilized, penalizing the number of active terms in the dynamic model following
the equation

$\displaystyle\operatorname*{arg\,min}_{\Xi}||\dot{\boldsymbol{x}}-\boldsymbol{\zeta}(\boldsymbol{x})\boldsymbol{\Xi}||_{2}+\boldsymbol{\nu}\boldsymbol{||\Xi||_{1}}$
(11)

with $\dot{\boldsymbol{x}}$ being the derivatives of the states and
$\boldsymbol{\zeta}(\boldsymbol{x})\boldsymbol{\Xi}$ a representation of the
dynamic system with candidate nonlinear functions
$\boldsymbol{\zeta}(\boldsymbol{x})$ and $\boldsymbol{\Xi}$ a coefficient
vector of the active terms. The use of $\boldsymbol{\Xi}$ in the regularizing
term $\boldsymbol{\nu}\boldsymbol{||\Xi||_{1}}$ accounts for a sparse model.
However, when trying to determine the unknown parameters
$\boldsymbol{\lambda}$ of a known set of equations, i.e. known set of
$\boldsymbol{\zeta}(\boldsymbol{x})$, the regularizing
$\boldsymbol{\nu}\boldsymbol{||\Xi||_{1}}$ can be neglected. For more
information, the interested reader can refer to [4]. In this paper, the
PySINDy library was used [12]. It should be mentioned here that SINDy is an
active research field, so improved approaches are proposed continuously; in
this paper, we went with the standard implementation in the PySINDy package.

##  III Results

To explore the performance of BPINNs, we assess them against two parameters in
this paper: first, for varying levels of measurement noise in the training
data, and, second, for varying lengths of the time series they receive as
training input. The noise level in data $\mathcal{D}$ (measurement noise)
follows an additive Gaussian noise model:

$\displaystyle\eta=\mathcal{N}(0,\gamma_{noise})$ (12)

with $\mathcal{N}$ being the Gaussian distribution of the measurement noise
and $\gamma_{noise}=\bar{\boldsymbol{|x|}}\cdot K$, with $K\in[0.0,0.05]$ and
$\bar{\boldsymbol{|x|}}$ the mean absolute value of the state.This means that
we examine noise levels that vary between 0% and 5%. As a second parameter, we
vary the length $T$ of the trajectories we use for training within the range
$T\in[12,27]\,\text{s}$ and we investigate its influence on the estimation
accuracy $\epsilon$.

For both cases, we measure the performance of BPINNs against SINDy and PINNs
using the percentage error (PE):

$\displaystyle\boldsymbol{\epsilon}=\left|\frac{(\hat{\boldsymbol{\lambda}}-\boldsymbol{\lambda})}{\boldsymbol{\lambda}}\right|\cdot
100\,\%$ (13)

with $\hat{\boldsymbol{\lambda}}$ representing the mean of the corresponding
posterior distribution (after 10 runs) in the case of BPINNs, and the mean of
the point estimate (again, after 10 runs) in the case of SINDy and PINNs. For
BPINNs, we also report the standard deviation
$\boldsymbol{\Gamma}_{\text{post}}$ of the posterior distributions for the
parameter estimates $\boldsymbol{\lambda}$, since BPINNs are able to deliver
the standard deviation of their expected estimate in contrast to SINDy and
PINNs. This can be interpreted as a measure of confidence for the estimates.
In Table II and Table III we report the normalized values of
$\boldsymbol{\Gamma}_{\text{post}}$, using the following equation:

$\displaystyle\boldsymbol{\tau}=\boldsymbol{\Gamma}_{\text{post}}\cdot\frac{100\%}{\boldsymbol{\lambda}}$
(14)

The most certain estimate would give a standard deviation
$\Gamma_{\text{post}}\approx 0$, and, therefore, $\tau\approx 0$. We expect
$\epsilon$ and $\Gamma_{\text{post}}$ to increase with increasing noise in
$\mathcal{D}$. At the same time, it is interesting to highlight that larger
deviations in the training trajectories $\boldsymbol{x}$ should lead to
smaller estimation errors $\epsilon$ and larger confidence (narrower
$\Gamma_{\text{post}}$), as there are less functions found by the BPINN that
fit $\mathcal{D}$. We will see the influence of the amount of deviations when
we examine the trajectory length in Section III-B

###  III-A Influence of noise

To investigate the performance of BPINN under noise, a noise level ranging
from $0\,\%$ to $5\,\%$ is added to each measurement, in steps of $1\,\%$
following (12). Table II compares the resulting performance of the BPINN to
the PINN and SINDy algorithm. We observe that all three algorithms deliver a
much more accurate estimate for the damping coefficient $d$ compared to the
inertia constant $m$. Let us highlight here that BPINNs offer richer
information than SINDy and PINNs through the standard deviation of the
posterior distribution $\Gamma_{\text{post}}$ (presented in its normalized
form $\tau$, see (14)). Through that, BPINNs allow us to understand that there
is a much higher degree of certainty in the estimate of $d$, where $\tau$ is
significantly smaller, compared to the estimate of $m$. The higher accuracy of
the $d$ estimate relates to $d$ being primarily dependent on the frequency
deviation $\Delta\omega$, while $m$ primarily depends on the rate of change of
frequency $\dot{\Delta\omega}$ (ROCOF). Inspecting Fig. 2, we observe that
frequency deviation $\Delta\omega$ varies more than ROCOF, and, thus, the
dataset contains more information useful for determining $d$. We continue this
analysis in the discussion of the results, see Section III-C.

Comparing the performance of the three methods, we observe that for $d$, BPINN
is not only the only one that offers a confidence measure, but also
outperforms both SINDy and PINNs in the point estimate, achieving the lowest
error $\epsilon_{d}$ across all noise levels and all fast/slow dynamics
scenarios. This is not the case, however, for the estimate of $m$. In that
case, PINNs appear to perform best for the fast dynamics, while SINDy performs
best for the slow dynamics, as long as there is no noise.

Besides the confidence measure that BPINNs offer, the most important result of
this investigation is the robustness of both BPINNs and PINNs to noise. With
noise gradually increasing from 0% to 5%, SINDy performance deteriorates
rapidly, reaching even a 50% estimation error. In contrast, BPINNs maintain an
almost constant very low error for the estimate of $d$, independent of the
noise level, while they also achieve a much lower error $\epsilon_{m}$ than
SINDy for higher noise levels. PINNs seem to be equally robust to noise as
BPINNs; and although they have higher errors than BPINNs for the estimate of
$d$, they do achieve the lowest error across all scenarios for the estimate of
$m$, with increasing noise levels.

TABLE II: Percent error (PE) $\epsilon$ and relative standard deviation $\tau$
for different noise levels

$K$ | Algor. | Fast dynamics 1 | Fast dynamics 2 | Slow dynamics 1 | Slow dynamics 2
---|---|---|---|---|---
|  | Inertia Error | Damping Error | Inertia Error | Damping Error | Inertia Error | Damping Error | Inertia Error | Damping Error
|  |  $\epsilon_{m}$ [$\%$] |  $\tau_{m}$ [$\%$] |  $\epsilon_{d}$ [$\%$] |  $\tau_{d}$ [$\%$] |  $\epsilon_{m}$ [$\%$] |  $\tau_{m}$ [$\%$] |  $\epsilon_{d}$ [$\%$] |  $\tau_{d}$ [$\%$] |  $\epsilon_{m}$ [$\%$] |  $\tau_{m}$ [$\%$] |  $\epsilon_{d}$ [$\%$] |  $\tau_{d}$ [$\%$] |  $\epsilon_{m}$ [$\%$] |  $\tau_{m}$ [$\%$] |  $\epsilon_{d}$ [$\%$] |  $\tau_{d}$ [$\%$]
$0\%$ | BPINN | 1.20 | 9.26 | 0.01 | 0.011 | 2.45 | 28.03 | 0.01 | 0.005 | 2.38 | 53.55 | 0.002 | 0.001 | 5.69 | 53.54 | 0.001 | 0.002
SINDy | 3.80 |  | 0.14 |  | 0.08 |  | 0.54 |  | 0.05 |  | 0.58 |  | 0.35 |  | 0.61 |
PINN | 0.34 |  | 0.84 |  | 0.19 |  | 1.22 |  | 1.08 |  | 1.41 |  | 3.16 |  | 1.15 |
$1\%$ | BPINN | 1.54 | 10.79 | 0.01 | 0.010 | 2.22 | 29.43 | 0.01 | 0.004 | 2.20 | 57.32 | 0.002 | 0.001 | 7.33 | 56.89 | 0.001 | 0.002
SINDy | 14.76 |  | 1.38 |  | 22.82 |  | 0.75 |  | 4.49 |  | 8.37 |  | 5.72 |  | 5.41 |
PINN | 2.23 |  | 1.98 |  | 1.62 |  | 1.44 |  | 1.65 |  | 1.60 |  | 2.53 |  | 1.23 |
$2\%$ | BPINN | 1.50 | 11.78 | 0.01 | 0.011 | 2.06 | 38.47 | 0.01 | 0.005 | 3.53 | 61.13 | 0.002 | 0.001 | 10.53 | 64.14 | 0.013 | 0.002
SINDy | 13.21 |  | 4.13 |  | 44.43 |  | 1.45 |  | 28.93 |  | 24.89 |  | 9.16 |  | 21.36 |
PINN | 0.11 |  | 2.46 |  | 0.91 |  | 2.01 |  | 1.73 |  | 1.74 |  | 3.85 |  | 1.22 |
$3\%$ | BPINN | 1.83 | 25.82 | 0.02 | 0.010 | 3.54 | 47.62 | 0.01 | 0.005 | 2.94 | 81.93 | 0.003 | 0.001 | 11.44 | 67.55 | 0.016 | 0.004
SINDy | 19.54 |  | 2.18 |  | 21.40 |  | 3.07 |  | 24.86 |  | 34.06 |  | 20.86 |  | 32.71 |
PINN | 2.71 |  | 2.90 |  | 1.23 |  | 2.35 |  | 2.77 |  | 1.77 |  | 3.69 |  | 1.48 |
$4\%$ | BPINN | 1.41 | 25.61 | 0.02 | 0.011 | 3.87 | 62.82 | 0.01 | 0.004 | 5.43 | 43.79 | 0.004 | 0.001 | 13.59 | 73.86 | 0.021 | 0.004
SINDy | 27.21 |  | 5.55 |  | 36.33 |  | 2.71 |  | 36.26 |  | 42.21 |  | 19.99 |  | 34.51 |
PINN | 0.77 |  | 3.04 |  | 0.61 |  | 2.95 |  | 2.89 |  | 1.83 |  | 4.36 |  | 1.42 |
$5\%$ | BPINN | 2.03 | 21.59 | 0.02 | 0.011 | 3.72 | 59.10 | 0.01 | 0.005 | 5.91 | 87.96 | 0.005 | 0.001 | 13.00 | 78.03 | 0.020 | 0.004
SINDy | 37.88 |  | 9.67 |  | 58.79 |  | 2.94 |  | 33.09 |  | 51.91 |  | 23.51 |  | 47.99 |
PINN | 0.41 |  | 3.99 |  | 2.29 |  | 3.26 |  | 3.27 |  | 2.06 |  | 2.07 |  | 1.52 |

###  III-B Influence of trajectory length

We expand our analysis, by including as parameters both the noise level and
the length of the time series that the three methods receive as training
input. Table III presents the results for varying levels of noise $K$ between
0% and 5%, similar to Section III-A, and varying trajectory length $T$ from
the default $27\,\text{s}$ down to $12\,\text{s}$.

Here again, we observe the excellent performance of BPINNs for the estimate of
$d$ across all noise levels and trajectory lengths, and the very good
performance of PINNs for the estimate of $m$, for higher noise levels and all
trajectory lengths. Here, it is interesting to highlight three phenomena.
First, that SINDy appears to perform better for shorter trajectory lengths
compared to the longer trajectory lengths in Table II, although still an order
of magnitude worse for higher noise, especially in systems with fast dynamics.
Second, PINNs appear to have a consistent good performance with errors up to
about 3% across all scenarios. Third, BPINNs appear to perform significantly
better in systems with faster dynamics compared to systems with slow dynamics.
We will examine the reasons behind these phenomena in Section III-C.

TABLE III: Percent error (PE) $\epsilon$ and relative standard deviation
$\tau$ for different noise levels and trajectory length

| Algor. | Fast dynamics 1 | Fast dynamics 2 | Slow dynamics 1 | Slow dynamics 2
---|---|---|---|---|---
($K$; |  | Inertia Error | Damping Error | Inertia Error | Damping Error | Inertia Error | Damping Error | Inertia Error | Damping Error
$T$) |  |  $\epsilon_{m}$ [$\%$] |  $\tau_{m}$ [$\%$] |  $\epsilon_{d}$ [$\%$] |  $\tau_{d}$ [$\%$] |  $\epsilon_{m}$ [$\%$] |  $\tau_{m}$ [$\%$] |  $\epsilon_{d}$ [$\%$] |  $\tau_{d}$ [$\%$] |  $\epsilon_{m}$ [$\%$] |  $\tau_{m}$ [$\%$] |  $\epsilon_{d}$ [$\%$] |  $\tau_{d}$ [$\%$] |  $\epsilon_{m}$ [$\%$] |  $\tau_{m}$ [$\%$] |  $\epsilon_{d}$ [$\%$] |  $\tau_{d}$ [$\%$]
$(0\%;$ | BPINN | 1.20 | 9.26 | 0.01 | 0.011 | 2.45 | 28.03 | 0.01 | 0.005 | 2.38 | 53.55 | 0.002 | 0.001 | 5.69 | 53.54 | 0.001 | 0.002
SINDy | 3.80 |  | 0.14 |  | 0.08 |  | 0.54 |  | 0.05 |  | 0.58 |  | 0.35 |  | 0.61 |
$27\,\text{s})$ | PINN | 0.34 |  | 0.84 |  | 0.19 |  | 1.22 |  | 1.08 |  | 1.41 |  | 3.16 |  | 1.15 |
$(1\%;$ | BPINN | 1.88 | 6.81 | 0.01 | 0.011 | 1.71 | 28.50 | 0.01 | 0.004 | 3.14 | 55.40 | 0.002 | 0.001 | 9.38 | 55.62 | 0.01 | 0.002
SINDy | 1.87 |  | 0.50 |  | 22.69 |  | 0.34 |  | 2.06 |  | 3.24 |  | 1.43 |  | 2.37 |
$24\,\text{s})$ | PINN | 1.13 |  | 1.29 |  | 0.53 |  | 0.535 |  | 2.01 |  | 1.56 |  | 2.01 |  | 1.28 |
$(2\%;$ | BPINN | 1.19 | 4.99 | 0.02 | 0.011 | 3.28 | 41.99 | 0.02 | 0.005 | 4.51 | 61.88 | 0.004 | 0.001 | 13.38 | 62.67 | 0.02 | 0.005
SINDy | 2.31 |  | 1.17 |  | 11.66 |  | 0.82 |  | 4.01 |  | 7.97 |  | 2.02 |  | 5.78 |
$21\,\text{s})$ | PINN | 1.96 |  | 1.27 |  | 0.32 |  | 1.58 |  | 0.92 |  | 1.52 |  | 1.08 |  | 1.24 |
$(3\%;$ | BPINN | 2.18 | 12.89 | 0.02 | 0.011 | 2.55 | 36.85 | 0.02 | 0.004 | 7.80 | 73.11 | 0.005 | 0.001 | 19.70 | 65.94 | 0.03 | 0.006
SINDy | 9.53 |  | 0.88 |  | 37.26 |  | 1.46 |  | 8.78 |  | 11.27 |  | 4.45 |  | 8.39 |
$18\,\text{s})$ | PINN | 0.23 |  | 1.18 |  | 0.83 |  | 1.72 |  | 1.14 |  | 1.72 |  | 1.84 |  | 1.47 |
$(4\%;$ | BPINN | 2.16 | 6.24 | 0.03 | 0.011 | 1.23 | 31.41 | 0.02 | 0.005 | 9.38 | 81.09 | 0.05 | 0.006 | 19.57 | 78.59 | 0.02 | 0.017
SINDy | 27.09 |  | 3.41 |  | 31.19 |  | 1.34 |  | 11.50 |  | 12.41 |  | 9.54 |  | 8.56 |
$15\,\text{s})$ | PINN | 0.71 |  | 1.42 |  | 0.32 |  | 2.72 |  | 2.21 |  | 2.43 |  | 2.42 |  | 2.20 |
$(5\%;$ | BPINN | 2.17 | 2.27 | 0.04 | 0.009 | 2.25 | 42.89 | 0.02 | 0.008 | 6.84 | 105.98 | 0.01 | 0.036 | 16.57 | 108.09 | 0.01 | 0.092
SINDy | 41.30 |  | 12.56 |  | 42.89 |  | 2.09 |  | 10.84 |  | 10.05 |  | 8.87 |  | 6.17 |
$12\,\text{s})$ | PINN | 2.04 |  | 2.24 |  | 1.58 |  | 1.40 |  | 2.85 |  | 1.65 |  | 3.05 |  | 1.45 |

###  III-C Discussion

#### BPINN vs SINDy Performance under noisy data

From Table II and Table III it is clear that in the presence of noise, the
BPINN provides estimates that are significantly better than SINDy, and
comparable to PINN. The reason for this can be found in the differentiation
methods. SINDy is calculating the derivative for every point in $\mathcal{D}$
and tries to find suitable parameters for the resulting differential
equations. So, in high noise, the SINDy derivatives can have large errors
between the single points. This is different for the BPINN and PINN, because
they fit the whole trajectory $\hat{\boldsymbol{x}}$ first, and then calculate
the derivative of the resulting fit $\frac{d}{dt}\hat{\boldsymbol{x}}$.

#### BPINN Performance – $m$ vs $d$

BPINNs appear to achieve 2 to 3 orders of magnitude lower estimation errors
for $d$ compared to $m$. This is more pronounced when estimating $m$ and $d$
in slow dynamics, and when using shorter trajectory lengths. Indeed, in slow
dynamics, the trajectory of $\boldsymbol{x}$ contains significantly less
information about $m$. To investigate the underlying reasons, in Fig. 3 we
plot the reconstruction of the $\Delta\omega$ trajectory for all input
trajectory lengths $T$ of the Slow Dynamics 2 scenario, using the BPINN
estimates found in Table III. We observe that the resulting trajectory is not
sensitive to varying parameter $m$: despite considering BPINN estimates which
had a non-negligible estimation error $\epsilon_{m}$, all trajectories look
similar and closely overlap with the true trajectory $x$. This also shows that
systems with higher inertia and damping are less sensitive to larger
variations of the inertia constant $m$. This, however, also results to larger
estimation errors, as the BPINN can fit a wider range of possible functions
from the family of functions to the training data. Same goes for
$\boldsymbol{\Gamma}_{post}$, as a larger $\boldsymbol{\Gamma}_{post}$ does
not result in significantly increased loss during training, since these
functions do not contradict the target trajectory. This is also the reason
that PINN performs better in these cases: PINNs have to fit only a single
function and not a family of functions as the BPINN does. Similar is the case
when we consider shorter trajectory lengths for training. Assuming the same
sampling frequency for all lengths $T$, there are fewer points that contradict
the trajectory estimate of the BPINN. So, the loss during training is small
for a wider range of functions, and this makes them appear as a suitable fit
for the optimizer. This is the reason that for lower noise, both SINDy and
PINN achieve more accurate estimates of $m$. Still, we shall mention that for
cases with slower dynamics, larger estimation errors of $m$ similar to BPINNs
do not seem to significantly affect the expected dynamic performance, as shown
in Fig. 3.

Figure 3: Reconstruction of $\Delta\omega$ with BPINN $m$ estimates for the
Slow Dynamics 2 scenario from Table III

#### Confidence measure

We claim that the magnitude of the normalized standard deviation $\tau$ not
only gives a form parameter of the posterior distribution but it also contains
information about how ‘confident’ the BPINN is that the delivered point
estimate, i.e. mean of posterior distribution, is accurate. In Fig. 4(b) we
plot the normalized standard deviation $\tau$ of the posterior distributions
against the error of the mean estimate $\epsilon$ for $m$ and $d$, for the
four scenarios. Indeed, we observe a clear correlation, where the standard
deviation is monotonously increasing with increasing estimation error for both
parameters. Such a confidence measure can be very valuable for the operation
of safety critical systems, such as power systems, as the magnitude of the
posterior distributions standard deviation can give a useful indication about
how much we can trust the point estimate of the BPINN about e.g. the system
inertia. Methods such as conventional neural networks or SINDy do not offer
any information about how ‘trustworthy’ is their output. Future work shall
perform detailed correlation analysis utilizing a wide range of scenarios.

(a) Standard deviation of all $m$ estimates versus the error

(b) Standard deviation of all $d$ estimates versus the error

Figure 4: Standard deviation $\tau$ versus estimate errors $\epsilon$

#### Initialization of the priors

Although BPINNs do not require prior knowledge about the target parameters or
the noise, they do require the initialization of the physical priors
$\boldsymbol{\lambda}$. It appears that this initialization should be chosen
carefully to achieve small errors $\epsilon$ when the noise in $\mathcal{D}$
is unknown, which would be the case when receiving measurement data. So, a
general prior that works for all noise levels and all scenarios should be
determined. We found through parameter studies that a normal distribution
$P(\boldsymbol{\lambda})$ can be used as follows

$\displaystyle P(\boldsymbol{\lambda})$
$\displaystyle\sim\mathcal{N}(1.0,5.0/P_{prec})$ (15) $\displaystyle P_{prec}$
$\displaystyle\sim\Gamma(1.0,0.1)$ (16)

where $P(\boldsymbol{\lambda})$ depends on a precision Gamma distribution of
$P_{prec}$ with shape parameter $\alpha=1$ and rate parameter $\beta=0.1$. The
resulting probability density function, besides being as uninformative as
possible due to the $\Gamma$ dependency, also limits the samples to being
positive. This is valid, since negative values for
$\boldsymbol{\lambda}=\\{m,d\\}$ are physically impossible. It also provides a
prior, that avoids a biased system, while providing the expected shape
(Gaussian), so the estimation result does not rely on previously available
knowledge, which the prior is normally utilized for.

#### Runtime

The average runtime of SINDy was $0.0032\,\text{s}$ for one estimate. The
BPINN training took on average $80\,\text{s}$ while PINN training took
$18.53\,\text{s}$. All methods were trained on an Intel i7 CPU. The increased
computation time for the BPINN can be explained by the family of functions
that has to be fit in comparison to the PINN, which needs to fit a neural
network to only a single function. Still, there is significant potential to
decrease the computation time for the BPINN with the utilization of GPU
computing and code optimization, which is object of our future work.

##  IV Conclusion and Outlook

This paper, to the best of our knowledge, introduced for the first time the
Bayesian Physics-Informed Neural Networks for applications in power systems.
In contrast to existing approaches, Bayesian Physics-Informed Neural Networks
(BPINNs) deliver a confidence measure of their estimates, while they do not
require inaccurate assumptions about the prior distributions of the variables
to be determined.

This paper applies the BPINNs for robust identification of the system inertia
and damping, using a single machine infinite bus system as the guiding
example. Comparing BPINNs with the Physics-Informed Neural Networks (PINNs)
and the recently popular approach for System Identification, SINDy, we assess
their performance. We find that BPINNs and PINNs are robust against all noise
levels, delivering estimates of the system inertia and damping with
significantly lower error compared to SINDy, especially as the noise levels
increased.

BPINNs combine the strengths of the Physics-Informed Neural Networks with a
confidence measure for their outputs, delivering a tool that can offer a
degree of ‘trustworthiness’ of their output. This opens up a wide range of
opportunities for applications in power systems, that remain to be explored.

## References

  * [1] D. Yang, B. Wang, G. Cai, Z. Chen, J. Ma, Z. Sun, and L. Wang, “Data-driven estimation of inertia for multiarea interconnected power systems using dynamic mode decomposition,” IEEE Transactions on Industrial Informatics, vol. 17, no. 4, pp. 2686–2695, 2021.
  * [2] Y. Susuki, R. Hamasaki, and A. Ishigame, “Estimation of power system inertia using nonlinear koopman modes,” in 2018 IEEE Power $\&$ Energy Society General Meeting (PESGM), pp. 1–5, 2018.
  * [3] J. Zhao, A. Gómez-Expósito, M. Netto, L. Mili, A. Abur, V. Terzija, I. Kamwa, B. Pal, A. K. Singh, J. Qi, Z. Huang, and A. P. S. Meliopoulos, “Power system dynamic state estimation: Motivations, definitions, methodologies, and future work,” IEEE Transactions on Power Systems, vol. 34, no. 4, pp. 3188–3198, 2019.
  * [4] S. L. Brunton, J. L. Procter, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” Applied Mathematics, vol. 113, no. 15, pp. 3932–3937, 2016.
  * [5] J. Stiasny, G. S. Misyris, and S. Chatzivasileiadis, “Physics-informed neural networks for non-linear system identification for power system dynamics,” in 2021 IEEE Madrid PowerTech, pp. 1–6, 2021.
  * [6] N. Petra, C. G. Petra, Z. Zhang, E. M. Constantinescu, and M. Anitescu, “A bayesian approach for parameter estimation with uncertainty for dynamic power systems,” IEEE Transactions on Power Systems, vol. 32, no. 4, pp. 2735–2743, 2017.
  * [7] L. Yang, X. Meng, and G. E. Karniadakis, “B-pinns: Bayesian physics-informed neural networks for forward and inverse pde problems with noisy data,” Journal of Computational Physics, vol. 425, p. 109913, 2021.
  * [8] M. Raissi, P. Perdikaris, and G. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” Journal of Computational Physics, vol. 378, pp. 686–707, 2019.
  * [9] I. Kononenko, “Bayesian neural networks,” Biological Cybernetics, vol. 61, no. 5, pp. 361–370, 1989.
  * [10] Q. Liu and D. Wang, “Stein variational gradient descent: A general purpose bayesian inference algorithm,” in Advances in Neural Information Processing Systems (D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett, eds.), vol. 29, Curran Associates, Inc., 2016.
  * [11] S. Stock, J. Stiasny, D. Babazadeh, C. Becker, and S. Chatzivasileiadis, “Online appendix: BPINN for system identification.” https://github.com/SimonStock/BPINN$\\_$for$\\_$System$\\_$Identification/, 2022\.
  * [12] A. Kaptanoglu, B. de Silva, U. Fasel, K. Kaheman, A. Goldschmidt, J. Callaham, C. Delahunt, Z. Nicolaou, K. Champion, J.-C. Loiseau, J. Kutz, and S. Brunton, “PySINDy: A comprehensive python package for robust sparse system identification,” Journal of Open Source Software, vol. 7, p. 3994, jan 2022.


# RME-GAN: A Learning Framework for Radio Map Estimation based on Conditional
Generative Adversarial Network

Songyang Zhang, Achintha Wijesinghe, and Zhi Ding * S. Zhang, and A.
Wijesinghe have the same contribution to this work.This work was supported in
part by the National Science Foundation under Grant 2029848 and Grant
2029027.S. Zhang, A. Wijesinghe and Z. Ding are with Department of Electrical
and Computer Engineering, University of California, Davis, CA, 95616. (E-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>and zding@ucdavis.edu).

###### Abstract

Outdoor radio map estimation is an important tool for network planning and
resource management in modern Internet of Things (IoT) and cellular systems.
Radio map describes spatial signal strength distribution and provides network
coverage information. A practical goal is to estimate fine-resolution radio
maps from sparse radio strength measurements. However, non-uniformly
positioned measurements and access obstacles can make it difficult for
accurate radio map estimation (RME) and spectrum planning in many outdoor
environments. In this work, we develop a two-phase learning framework for
radio map estimation by integrating radio propagation model and designing a
conditional generative adversarial network (cGAN). We first explore global
information to extract the radio propagation patterns. We then focus on the
local features to estimate the effect of shadowing on radio maps in order to
train and optimize the cGAN. Our experimental results demonstrate the efficacy
of the proposed framework for radio map estimation based on generative models
from sparse observations in outdoor scenarios.

###### Index Terms:

Radio map estimation (RME), conditional generative adversarial networks
(cGAN), radio measurement, network planning.

##  I Introduction

The rapid development of information technologies has spawned many novel
concepts and applications, such as the Internet of Things (IoT), autonomous
driving, and cellular networks. For these wireless services, efficient
spectrum planning and network coverage analysis are critical to spectrum
efficiency and user experience. Specifically, radio maps provide vital
information on spatial distribution of radio frequency (RF) signal power, and
facilitate resource allocation and network outage diagnosis [1]. In general,
RF ‘radio map’ refers to a geographical signal power spectrum density (PSD),
formed by cumulative received RF signal power, as a function of spatial
locations and frequencies, and may also exhibit temporal variations [2].
Enriching useful information on propagation behavior and spectrum occupancy,
radio map could provide more detailed information on PSD distribution [3]
spatially, shown as Fig. 1(a). However, RF signal powers can only be sparsely
measured in practice. Thus, a practical challenge is the estimation of fine-
resolution radio map from spatially sparse signal power samples collected by
sensors or user devices as shown in Fig. 1(b). Our goal here is to achieve
more efficient and accurate estimation of radio maps from sparse observations.

(a) Radio Map over Different Locations

(b) Sparse Samples Collected by Sensors

Figure 1: Example of Radio Map: a) radio map generated from simulators with
three transmitters; b) 1000 samples from the original radio map.

Most current approaches to radiomap estimation (reconstruction) can be
categorized into either model-based or model-free. Model-based methods usually
assume certain signal propagation models and express the received signal PSD
as a combination from active transmitters [2]. For example, the authors of [4]
introduced log-distance path loss model (LDPL) for Wi-Fi radio map
reconstruction. In [5], another model-based method proposed interpolation
based on thin-plate splines kernels. On the other hand, model-free methods do
not assume a specific model but explore more neighborhood information. Typical
examples include Radial Basis Function (RBF) [6] interpolation and inverse
distance weighted (IDW) interpolation [7]. In addition to interpolation
methods, there have also been a number of works that apply machine learning
for radio map reconstruction. The authors of [8] presented a fast radio map
estimation based on Unet [9]. Another team [10] proposed deep auto-encoders
for estimating geographical PSD distribution. We will provide a more complete
literature review later in Section II.

Despite many successes, existing solutions still exhibit certain limitations.
First, unlike the problem of astronomical image analysis [11] and 3D city map
reconstruction [12], ground truth radio map is more sensitive to the complex
environment beyond mere radio propagation models. Failure to capture effects
from shadowing or obstacles limits the efficacy of model-based approaches. On
the other hand, the performance of model-free methods relies on the quality of
observed samples. Generally, they assume a uniform distribution of
measurements over the geometric space. However, in practical applications such
as minimization of drive test (MDT) [13], measurements are collected from user
devices, which are typically non-uniformly distributed within cell coverage.
It is a challenge to process unevenly distributed observations efficiently.
Moreover, samples in different training sets do not necessarily have the same
radio propagation models and parameters. In view of their respective strengths
and shortcomings, we should consider effective integration of model-based and
model-free methods.

To address the aforementioned problems, we investigate radio propagation in
the outdoor scenario, and propose a two-phase learning framework (RME-GAN) for
radio map estimation based on conditional generative adversarial networks
(cGAN) [14]. Since a radio map can be dynamically constructed from real-time
measures and continuously updated [2], we mainly focus on the step of radio
map estimation from sparse observations over geographical locations. In phase
1, we integrate learning-based and model-based principles to extract global
information while leveraging radio propagation model. In phase 2, we utilize
different sampling strategies to capture local radiomap features and to
explore the shadowing effect. More specifically, we propose to capture the
sharp features in the Fourier domain for detailed information, before
resampling the measured signals geometrically to address the non-uniformity
among distributed observations. We summarize our contributions as follows:

  * •

We introduce a strategy to integrate well known radio propagation model into
data-driven approaches to guide the learning machines for global radiomap
information in the first phase of our proposed framework;

  * •

We design different sampling strategies to capture local radio information and
shadowing effects to overcome the non-uniformity of distributed observations
in the second phase of the proposed framework;

  * •

We provide a comprehensive experimental comparison with state-of-the-art
methods in terms of both reconstruction accuracy and fault diagnosis;

The experimental results demonstrate the ability of the proposed strategies to
capture global model information and local shadowing features, as well as the
power of the two-phase learning frameworks in radio map reconstruction.

We organize the remainder of this work as follows. We first overview the
related works of radio map reconstruction and generative adversarial network
(GAN) models in Section II. We then introduce our problem formulation in
Section III. In Section IV, we present the details of the proposed two-phase
learning frameworks based on cGAN. Following the presentation of the
experimental results in Section V, we conclude our works in Section VI.

##  II Related Work

Here we provide an overview of radio map estimation from sparse samples based
on either model-based or model-free approaches. We also introduce GAN models.

###  II-A Radio Map Estimation

####  II-A1 Model-based Radio Map Estimation

Model-based approaches usually leverage a certain signal propagation model.
Specifically, a radio map function is modeled as a function of frequency $f$
and a geographic location $\mathbf{c}$ by

$\mathbf{\Phi}(\mathbf{c},f)=\sum_{i=1}^{N_{t}}g_{i}(\mathbf{c},f)\Psi_{i}(f)$
(1)

where $N_{t}$ denotes the number of transmitters, $g_{i}$ denotes the channel
power gain of the $i$th transmitter, and $\Psi_{i}(f)$ is the PSD of the the
$i$th transmission [2]. Here, $g_{i}$ can either be assumed as an explicit
function with prior knowledge depending on specific tasks, or be represented
in terms of sensor parameters. For example, the authors of [4] modeled a
single narrowband WiFi transmitter with $g_{i}$ according to the log-distance
path loss model (LDPL). In [5], the authors modeled $g_{i}$ as a kernel
expansion to estimate the power gain. Other model-based methods include
parallel factor analysis (PARAFAC) [15] and fixed rank kriging [16].

####  II-A2 Model-free Radio Map Reconstruction

Different from model-based methods, model-free methods usually do not assume a
specific propagation model but favor neighborhood information. In general, the
model-free methods can be categorized as interpolation-based methods and
learning-based methods as follows:

  * •

Interpolation-based Methods: Interpolation-based approaches express the PSD at
a particular location as a combination of neighborhood measurements, i.e.,

$\mathbf{\Phi}(\mathbf{c},f)=\sum_{i=1}^{N_{s}}w_{i}(\mathbf{c},f)q_{i}(f),$
(2)

where $q_{i}$ denotes the observation from the $i$th observation and $w_{i}$
is the combination weight [2]. One typical example is inverse distance
weighted (IDW) interpolation [7], where the power gain is inversely
proportional to the distance between transmitters and receivers. Besides
linear interpolation, an effective alternative is radio basis function (RBF)
interpolation [6], where different kernels, such as Gaussian, multiquadrics,
or spline, can be applied. Another classic interpolation approach is ordinary
Kriging methods [17], where the PSD is optimized by minimizing the variance of
estimation error. Beyond these traditional methods, recent works of radio map
interpolation studied integration with novel data analysis concepts, such as
graph signal processing [18], crowdsourcing methods [19], image inpainting
[20, 21] and tensor completion [22].

  * •

Learning-based Methods: Recent applications and success of machine learning in
wireless networks and IoT intelligence have stimulated more interests in
learning-based radio map estimation as another promising direction. Different
from explicit interpolation, learning-based methods focus on finding a direct
mapping from an input geometric map to its output PSD measurement. The
functional relationship between input and output is usually captured by a
neural network. For example, Unet is introduced in [8] for predicting the
pathloss. Similarly, a deep auto-encoder is another design for the deep neural
network mapping [10]. Other learning frameworks for radio map estimation
include transfer learning [23], generative adversarial network [24, 25, 26],
deep Gaussian process [27], reinforcement learning [28], and deep neural
networks [29, 30, 31].

Interested readers could refer to [32] for a more thorough overview of data-
driven radio map reconstruction.

###  II-B Generative Adversarial Networks

The idea of Generative Adversarial Networks (GAN) is to study a collection of
training examples to acquire sufficient information on sample distribution in
order to generate data samples with similar distribution [33]. In GAN, two
neural networks interact in the form of a zero-sum game, in which the
generator $G$ aims to generate new data with the same statistics as the
training dataset while the discriminator $D$ focuses on the identification of
the true and fake data samples. Suppose the generator learns the distribution
$p_{g}$ over data $\bf{x}$ by mapping a prior noise distribution
$p_{z}(\bf{z})$ to data space. The GAN model can be trained via

$\displaystyle\min_{G}\max_{D}\quad V(D,G)=\mathbb{E}_{x\sim p_{\rm
data}(\bf{x})}[\log D(\bf{x})]+$ $\displaystyle\mathbb{E}_{z\sim
p_{z}(\bf{z})}[\log(1-D(G(\bf{z})))],$ (3)

where $\mathbb{E}[\cdot]$ refers to the expectation operation.

Figure 2: Example of a Customized cGAN.

With the development of machine learning, GAN can be extended to a conditional
model (cGAN) if the generator and discriminator are conditioned on some extra
information $\mathbf{I}$ [14]. Then, the objective function could be set as

$\displaystyle\min_{G}\max_{D}\quad V(D,G)=\mathbb{E}_{x\sim
p_{data}(\bf{x})}[\log D(\bf{x|I})]+$ $\displaystyle\mathbb{E}_{z\sim
p_{z}(\bf{z})}[\log(1-D(G(\bf{z|I})))].$ (4)

Note that, previous works have shown that cGAN can be customized to learn the
desired properties by mixing the traditional loss function with cGAN
objectives [34] and utilizing the corresponding feature vectors as input. See
the illustration in Fig. 2. For example, to generate a painted image from
sketch, one can input the sketch to the generator and formulate the training
process as a one-to-one mapping. No random noise vector should be involved for
such mapping problem. The customized cGAN frameworks have achieved successes
in applications, including image painting [35], image-to-image translation
[36], and image manipulation [37].

##  III Problem Formulation

Now, we present the problem formulation in this work. In general, we focus on
the radio map reconstruction from sparse observations in the outdoor scenario.

Our model considers a set of regions with sparse observations of network
coverage, similar to Fig. 1(b). Each region is conformed to a size of
$256\times 256$ grid and contains $N_{t}$ transmitters with known positions.
Let $\mathbf{P_{t}}_{i}\in\mathbb{R}^{Nt\times 2}$ denote the 2-dimensional
(2D) coordinates of transmitters for the $i$th region. Region $i$ is also
attributed with their corresponding urban maps
$\mathbf{m}_{i}\in\mathbb{R}^{256\times 256}$ with building information.
Instead of knowing the fine-resolution radio map of the whole space, each
region only observes a set of $K_{i}$ sparse samples
$(\mathbf{x}_{i},\mathbf{c}_{i})$ collected from sensors, where $K_{i}\leq
256\times 256$, $\mathbf{x}_{i}\in\mathbb{R}_{+}^{K_{i}}$ are the non-negative
observed PSD, and $\mathbf{c}_{i}\in\mathbb{R}^{K_{i}\times 2}$ are the 2D
coordinates of the corresponding observations. Thus, each region $i$ is
characterized by

$\mathcal{F}_{i}=\\{(\mathbf{x}_{i},\mathbf{c}_{i}),\mathbf{m}_{i},\mathbf{P_{t}}_{i}\\}.$
(5)

Suppose that some regions have the ground truth of the fine-resolution radio
map $\mathbf{y}$. Treating these regions as training regions, our goal is to
train a learning machine $g(\cdot)$ to estimate the fine-resolution radio map
for all regions, i.e.,

$\tilde{\mathbf{y}}_{i}=g(\mathcal{F}_{i})\in\mathbb{R}_{+}^{256\times 256}.$
(6)

Unlike existing learning-based radio map reconstruction, we do not apply
strict constraints on the distribution of samples and propagation model
parameters. The following conditions are relaxed:

  * •

The sparse observations $\mathbf{x}_{i}$ can be either uniformly or unevenly
distributed over the whole region;

  * •

Different regions may have different number $K_{i}$ of observations;

  * •

The transmitters in each region may have different radio propagation
parameters.

To deal with such non-uniformly distributed observations and the unbalanced
data samples over different training regions, we propose a two-phase learning
framework for radiomap estimation (RME) based on GAN, known as RME-GAN. Our
goal is to train an efficient generator $g(\cdot)$ to estimate the full radio
map as described in Eq. (6).

##  IV Radio Map Reconstruction via RME-GAN

In this section, we present the design of the proposed RME-GAN, with the
overall structure shown as Fig. 3. We first introduce the basic sketch of RME-
GAN in Section IV-A, before describing the two-phase customization of features
and loss functions in Section IV-B and Section IV-C, respectively.

(a) Structure of Generator

(b) Structure of Discriminator

Figure 3: Structure of RME-GAN: a) The generator takes input $\mathcal{F}_{i}$ and consider the model information during training, where two phases of loss functions are applied to train the generator; b) The discriminator consists of 4 layers and generates a one-dimensional probability output. The discriminator input consists of 3 channels: the generated/true radio map (marked as orange), and the prior information as two-channel one-hot encodings (marked as grey and black).  TABLE I: Structure of Unet Layer | Input | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
---|---|---|---|---|---|---|---|---|---|---
Resolution | 256 | 256 | 128 | 64 | 64 | 32 | 32 | 16 | 8 | 4
Channel | 3 | 6 | 40 | 50 | 60 | 100 | 100 | 150 | 300 | 500
Filter | 3 | 5 | 5 | 5 | 5 | 3 | 5 | 5 | 5 | 4
Layer | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | Output
Resolution | 8 | 16 | 32 | 32 | 64 | 64 | 128 | 256 | 256 | 256
Channel | 300+300 | 150+150 | 100+100 | 100+100 | 60+60 | 50+50 | 40+40 | 20+6+3 | 20+3 | 1
Filter | 4 | 4 | 3 | 6 | 5 | 6 | 6 | 5 | 5 | -

###  IV-A Sketch of RME-GAN

RME-GAN learns a mapping function $g(\cdot)$ from the region features
$\mathcal{F}_{i}$ to its output estimated full radio map
$\tilde{\mathbf{y}}_{i}$. Since it is also a one-to-one mapping problem, we
could design RME-GAN based on the general cGAN structure similar to Fig. 2.
More specifically, we utilize an 18-layer deep Unet [8] structure for the
sketch of the generator, with parameters shown in Table. I.

To reformat the input features in $\mathcal{F}_{i}$, we first interpolate the
sparse observations $(\mathbf{x}_{i},\mathbf{c}_{i})$ to the same size as
region size via zero padding, i.e.,
$\mathbf{x}_{i}^{\prime}\in\mathbb{R}_{+}^{256\times 256}$. Since we consider
a convolutional neural network (CNN) structure for Unet, the interpolated zero
would not influence other features. We then encode the transmitter locations
$\mathbf{P_{t}}_{i}$ also in a $256\times 256$ matrix $\mathbf{p}_{i}$, where
the transmitter positions in $\mathbf{p}_{i}$ are marked as $1$ and other
positions are $0$. Features $\mathcal{F}_{i}$ are encoded to three channels of
features with size $256\times 256$, i.e.,
$\mathbf{f}_{i}=\\{\mathbf{x}_{i}^{\prime},\mathbf{m}_{i},\mathbf{p}_{i}\\}$.

With the encoded features, the basic objective function of RME-GAN can be
designed by

$\displaystyle\min_{G}\max_{D}\quad V(D,G)=\mathbb{E}_{y,I\sim p_{\rm
data}(\bf{y},\bf{I})}[\log D(\bf{y,I})]+$ $\displaystyle\mathbb{E}_{f\sim
p_{\rm data}(\bf{f})}[\log(1-D(G(\bf{f}),\mathbf{I}))],$ (7)

where $p_{\rm data}(\bf{y,I})$ is the joint distribution with the prior
knowledge $\mathbf{I}$. More specifically, $\mathbf{I}$ can be designed as a
function of the region features, i.e., $\mathbf{I}=f(\mathcal{F}_{i})$ or some
determined labels, such as fake or real [35]. Here, we set
$\mathbf{I}\in\mathbb{R}^{256\times 256\times 2}$ to represent the one-hot
encodings in two-channel to give prior information to the discriminator as the
input is a generated or real radio map.

We can split the min-max optimization problem to sub-problems, and train the
generator $G$ and discriminator $D$ iteratively. Specifically, the loss
function to optimize $D$ is

$L_{D}=-V(D,G),$ (8)

and the general loss function to optimize $G$ is

$L_{G}=\mathbb{E}_{x\sim p_{\rm
data}(\bf{x})}[\log(1-D(G(\bf{x}),\mathbf{I}))].$ (9)

The overall sketch of RME-GAN is illusrated in Fig. 3. Note that the cGAN can
be customized by adding additional loss terms to capture desired properties.
In this work, we consider a two-phase loss function design for Eq. (9) to
capture the global and local information, respectively. More details will be
discussed in Section IV-B and IV-C.

###  IV-B Two-phase Customization in RME-GAN

Earlier, we have identified two challenges in radio map estimation: 1) how to
efficiently integrate radio propagation model with learning machines to
capture global information; and 2) how to estimate the details from biased and
non-uniformly distributed observations. The first challenge lies in the
estimation of global model information while the second challenge deals with
shadowing effect and obstacle impact. To handle these challenges, we utilize
different sampling strategies to extract global patterns and local features,
respectively, as follows.

####  IV-B1 Upsampling for global information

In general, radio propagation should follow certain pathloss models, such as
LDPL [4]. Thus, to model global pathloss, we first upsample the sparse
observations $(\mathbf{x}_{i},\mathbf{c}_{i})$ to a radio map template
$\mathbf{z}\in\mathbb{R}_{+}^{256\times 256}$ via model-based interpolation
(MBI). Suppose that a system has $N_{t}$ transmitter with positions
$(a_{k},b_{k})$, $1\leq k\leq N_{t}$. In MBI, the signal strength (in dB) at
position $(a,b)$ can be modeled as

$P(a,b)=\sum_{k=1}^{N_{t}}(\alpha_{k}-10\theta_{k}\log_{10}(\sqrt{(a_{k}-a)^{2}+(b_{k}-b)^{2}})),$
(10)

where $\alpha_{k}$ is the power from $k$th transmitter at the reference
distance. Then, we optimize the paramters of the LDPL model from the
observations by

$\min_{\alpha_{k},\theta_{k}}\sum_{j=1}^{K_{i}}(\mathbf{x}(j)-P(\mathbf{c}(j)))^{2}$
(11)

where $\mathbf{x}(j)$ is the $j$th observation with position $\mathbf{c}(j)$.
By apply the LDPL model with estimated parameters, we can upsample the sparse
observations $\mathbf{x}_{i}$ as $\mathbf{z}_{i}\in\mathbb{R}_{+}^{256\times
256}$.

Although MBI may be not able to capture full effect of obstacles and
shadowing, upsampled template can still capture the global propagation pattern
together with MBI model parameters. Consideration of upsampled features during
the training process can provide a general guideline on radio map
reconstruction. We introduce how to use $\mathbf{z}_{i}$ to design loss
functions of Phase 1 in Section IV-C.

####  IV-B2 Downsampling for detailed features

After estimating a template of reconstructed radio maps, we downsample the
sparse observations in the second phase to correct errors in the details. We
consider two downsampling strategies:

  * •

Geometry-wise downsampling: As discussed in Section III, sparse observations
may be non-uniformly distributed over the whole space. To deal with unevenly
distributed samples, we apply geometric information to downsample the
observations. To balance the sample numbers in different locations, we first
build superpixels [38] based on urban map or interpolated template to split
the entire region into $N_{s}$ pieces. In each piece, we only pick one
observation of peak energy. Then, we have a downsampled version of
observations $\mathbf{x}\in\mathbb{R}_{+}^{K}$ as
$\mathbf{x_{s}}\in\mathbb{R}_{+}^{N_{s}}$. In addition to superpixel-based
space cut, other methods such as Voronoi diagram [4] or spectral clustering
[39] can be also applied.

  * •

Frequency-wise downsampling: In digital signal processing, details and
outliers are related to high frequencies in Fourier space. Thus, we could also
transform the estimated radio map $\tilde{\mathbf{y}}$ to Fourier space as
$\hat{\mathbf{y}}$. Then, we can select the first $N_{f}$ Fourier coefficients
in the high frequency part as additional information in loss functions of the
second phase, i.e., $\hat{\mathbf{y}}_{f}\in\mathbb{R}^{N_{f}}$.

With the extracted features under different sampling strategies, we propose a
two-phase design on the loss functions for the training of the generator,
i.e.,

$L_{G}^{\prime}=\left\\{\begin{aligned}
&L_{G}+L_{Global}\quad&\text{if}\quad\epsilon<\tau\\\
&L_{G}+L_{Local}\quad&\text{if}\quad\epsilon\geq\tau\end{aligned}\right.,$
(12)

where $L_{Global}$ consist of loss functions favoring global model information
and $L_{Local}$ are the loss functions focusing on the local details. Here, a
decision variable $\epsilon$ can be designed according to certain criteria,
i.e., length of epochs or validation accuracy. For convenience, we keep a
validation data set to evaluate the validation accuracy at each epoch to
determine the phase change point. The specific design of loss functions will
be introduced next in Section IV-C.

###  IV-C Loss Functions

After introducing the overall frameworks, we now present the design of loss
functions.

####  IV-C1 Phase 1 - Global Model Prediction

We first introduce the design of loss function in the first phase.

  * •

Baseline loss function: The first term is the mean squared error (MSE) between
the generated radio map and the ground truth $\mathbf{y}$:

$L_{\rm MSE}=\frac{1}{256\times
256}\cdot||\mathbf{y}-\tilde{\mathbf{y}}||_{F}^{2},$ (13)

where $||\cdot||_{F}$ is the Frobenius norm and
$\tilde{\mathbf{y}}=g(\mathbf{f})$. We also consider a second term to be the
total variation which describes the smoothness of the generated radio map
[35]. More specifically, we apply the total variation of images

$L_{\rm
TV}=\sum_{i,j\in\mathcal{N}}||\tilde{\mathbf{y}}(i)-\tilde{\mathbf{y}}(j)||_{2}^{2},$
(14)

where $\mathcal{N}$ indicates pixel neighborhood and $\tilde{\mathbf{y}}(i)$
is the $i$th element of the estimated radio map $\tilde{\mathbf{y}}$.

  * •

Loss function of global information: In addition to the baseline loss
functions, we also utilize upsampled radio map $\mathbf{z}$ from MBI to help
the learning of global radio propagation patterns. Instead of directly
evaluating the difference between $\tilde{\mathbf{y}}$ and $\mathbf{z}$, we
focus on gradient patterns. Such gradient patterns feature smoothness of radio
propagation from transmitters to the surroundings. Since the first phase aims
to learn the radio propagation model, MBI can help RME-GAN to learn the
propagation model and speed up training. To obtain the gradient-based loss
function, we first calculate the gradient of each pixel in the radio map for
up, down, right and left, i.e.,

$\mathcal{G}(\tilde{\mathbf{y}}(i))=[\mathcal{G}_{\rm up},\mathcal{G}_{\rm
down},\mathcal{G}_{\rm right},\mathcal{G}_{\rm left}]\in\mathbb{R}^{4}.$ (15)

Then, we calculate the cosine similar between
$\mathcal{G}(\tilde{\mathbf{y}}(i))$ and $\mathcal{G}({\mathbf{z}}(i))$. Then,
the loss function of global information is defined as

$L_{\rm Gradient}=\sum_{i=1}^{256\times
256}CS(\mathcal{G}(\tilde{\mathbf{y}}(i)),\mathcal{G}({\mathbf{z}}(i))),$ (16)

where $CS(\cdot)$ refers to the cosine similarity [40].

####  IV-C2 Phase 2 - Local Shadowing Estimation

Next, we introduce the loss functions for local features in the second phase.

  * •

Baseline loss function: Similar to the first phase, the generated radio map
should be smooth over the entire space. Thus, we also use total variation
$L_{\rm TV}$ as one of the baseline loss functions. In addition, the $L_{\rm
MSE}$ is also used here to control the overall accuracy. However, to emphasize
more details depending on the downsampled observations, we assign a much
smaller weight on $L_{\rm MSE}$ in the second phase compared to the first
phase.

  * •

Geometric loss function: To mitigate the effect of non-uniformly distributed
observations, we favor downsampled observations among all observations. Since
the sparse observations are also true power gain in the radio map, we define a
geometric loss function as the MSE between downsampled
$\mathbf{x}_{s}=\mathbf{y}_{s}$ and the generated PSD at the same positions,
denoted by $\tilde{\mathbf{y}}_{s}$, i.e.,

$L_{\rm Geo}=MSE(\mathbf{x}_{s},\tilde{\mathbf{y}}_{s}).$ (17)

  * •

Frequency domain loss function: As discussed before in Section IV-B, high
frequency coefficients in the Fourier space can capture the details or
outliers in a radio map. To capture such information, we also calculate the
MSE between $\hat{\mathbf{y}}_{f}$ and $\hat{\tilde{\mathbf{y}}}_{f}$ denoted
by

$L_{\rm HPF}=MSE(\hat{\mathbf{y}}_{f},\hat{\tilde{\mathbf{y}}}_{f}),$ (18)

where $\hat{\mathbf{y}}_{f}$ includes the first top $N_{f}$ high frequency
coefficients in the Fourier domain.

  * •

Multi-scale structural similarity index: In image processing, multi-scale
structural similarity index (MS-SSIM) is a useful term to evaluate the image
quality based on human visual system [41]. Since details in radio map can be
viewed as a fine-resolution image while the model-based template can be viewed
as a coarse-resolution image, we could also consider the MS-SSIM index to
improve the quality of generated radio map in the second phase. Let
$\mathbf{u}={u_{i}|i=1,2,\cdots,N}$ and $\mathbf{v}={v_{i}|i=1,2,\cdots,N}$ be
two discrete non-negative signals, e.g., two image patches from the same
position or two radio map residing on the same region. The luminance, contrast
and structure comparison measurements can be calculated as follows:

$l(\mathbf{u},\mathbf{v})=\frac{2\mu_{u}\mu_{v}+C_{1}}{\mu_{u}^{2}+\mu_{v}^{2}+C_{1}},$
(19)
$c(\mathbf{u},\mathbf{v})=\frac{2\sigma_{u}\sigma_{v}+C_{2}}{\sigma_{u}^{2}+\sigma_{v}^{2}+C_{2}},$
(20)
$s(\mathbf{u},\mathbf{v})=\frac{2\sigma_{uv}+C_{3}}{\sigma_{u}\sigma_{v}+C_{3}},$
(21)

where $\mu$ represents mean, $\sigma$ is the variance, $C_{i}$ are small
constants introduced in [42]. To calculate MS-SSIM, the system iteratively
applies a low-pass filter and downsamples the filtered image by a factor of 2
on reference and distorted image signals as inputs. Suppose that $M$ is the
highest scale. The overall MS-SSIM is calculated by

${\rm
SSIM}(\mathbf{u},\mathbf{v})=[l_{M}(\mathbf{u},\mathbf{v})]^{\alpha_{M}}\cdot\prod_{j=1}^{M}[c_{j}(\mathbf{u},\mathbf{v})]^{\beta_{j}}[s_{j}(\mathbf{u},\mathbf{v})]^{\gamma_{j}}.$
(22)

Here, ${\rm SSIM}(\mathbf{u},\mathbf{v})\leq 1$, and ${\rm
SSIM}(\mathbf{u},\mathbf{v})=1$ iff $\mathbf{u}=\mathbf{v}$. The MS-SSIM loss
function is defined by

$L_{\rm SSIM}=1-{\rm SSIM}(\mathbf{y},\tilde{\mathbf{y}}).$ (23)

More properties of MS-SSIM can be found in [41].

###  IV-D Summary of Two-Phase Training Strategy

Combining all the loss terms aforementioned, the final loss function to train
the generator is represented as

$\displaystyle L_{G}^{\prime}=\left\\{\begin{aligned}
&\lambda_{1}L_{G}+\lambda_{2}L_{\rm MSE}+\lambda_{3}L_{\rm
TV}+\lambda_{4}L_{\rm Gradient},&\epsilon<\tau\\\
&\lambda_{1}L_{G}+\lambda_{2}L_{\rm MSE}+\lambda_{3}L_{\rm TV}+\\\
&\lambda_{5}L_{\rm SSIM}+\lambda_{6}L_{\rm Geo}+\lambda_{7}L_{\rm
HPF}&\epsilon\geq\tau\end{aligned}\right.$

where $\lambda_{i}$ are the weights to control each loss term. In the first
phase, we aims to apply model-based interpolation to help RME-GAN learn a
template to capture the radio propagation model. After enough epochs, we favor
more on the local information to correct errors in the details on the template
by learning the shadowing effect and obstacle impact. Although some baseline
loss functions are shared in both phases, the weights can be different
depending on the specific tasks.

##  V Experimental Results

In this section, we present the experimental results of the proposed methods.

###  V-A Experiment Setup

(a) Example of Urban Images.

(b) Example of Radio Map

Figure 4: Examples of RadioMapSeer Dataset: a) Urban maps with buildings as
white and background as black; b) radio map with pathloss as yellow.

(a) Ground Truth

(b) Three Selected Regions

Figure 5: One Example of Radio Map

####  V-A1 Dataset

In this work, we test the performance of the proposed method in the
RadioMapSeer dataset111https://radiomapseer.github.io/ [8]. The RadioMapSeer
dataset consists of 700 maps with 80 transmitter locations per map. The city
maps are taken from OpenStreetMap [43] covering metropolitans such as Ankara,
Berlin, Glasgow, Ljubljana, London and Tel Aviv. The heights of the
transmitters, receivers and buildings are set to 1.5m, 1.5m, and 25m,
respectively, in this dataset. The transmitter power equals 23dBM and the
carrier frequency is 5.9GHz. The radio maps are generated using the software
WinProp [44], and are saved in a 2D grid of $256\times 256m^{2}$. The radio
map has a 1-meter resolution. More specifically, we test on the higher
accuracy simulations in the RadioMapSeer datasets. An example of the
RadioMapSeer dataset is shown in Fig. 4, where building segmentation are
consider as the urban map in our framework. The sparse observations are
sampled from each radio map for training and testing.

(a) NMSE

(b) RMSE

Figure 6: Performance of Different Methods over Different Ratio of
Observations.

####  V-A2 Data Splitting

As mentioned in Section III, the sparse observations can be either uniformly
or non-uniformly distributed over the whole space, and different regions may
have different number of observations. Thus, we sample the sparse observations
under three different scenarios:

  * •

Setup 1 (Uniformly distributed): We uniformly sample $1\%$ of sparse
observations from the whole radio map.

  * •

Setup 2 (Unbalanced sample distribution among regions): In this setup, we
sample the radio map at a random ratio from $1\%\sim 10\%$ in each region, and
the sparse observations are uniformly distributed.

  * •

Setup 3: (Non-uniformly distributed) We divide the radio map into two sides
randomly. Then, we sample $1\%$ of observations from one side and $10\%$ from
the other side.

We compare the proposed RME-GAN with other methods in all the three scenarios
to demonstrate its efficacy. The whole dataset consists of 700 regions. For
each experiment scenario, we split the dataset with 500 regions as training
data, 100 regions for testing, and 100 regions for validation.

TABLE II: Performance in Different Setups | Setup 1 | Setup 2 | Setup 3
---|---|---|---
Model | NMSE | RMSE | NMSE | RMSE | NMSE | RMSE
Kriging | 0.8863 | 0.1947 | 0.8987 | 0.1962 | 0.9118 | 0.1971
RBF | 0.3532 | 0.1343 | 0.1830 | 0.0884 | 0.3215 | 0.1295
MBI | 0.8837 | 0.1973 | 0.5158 | 0.1499 | 0.5174 | 0.1500
cGAN | 0.1559 | 0.0904 | 0.15077 | 0.0888 | 0.1506 | 0.0889
AE | 0.1823 | 0.0978 | 0.2885 | 0.1238 | 0.2783 | 0.1217
Deep AE | 0.1898 | 0.0998 | 0.3152 | 0.1295 | 0.3058 | 0.1274
Unet | 0.0093 | 0.0220 | 0.0053 | 0.0166 | 0.0050 | 0.0161
RadioUnet | 0.0052 | 0.0164 | 0.0042 | 0.0148 | 0.0046 | 0.0155
RME-GAN | 0.0043 | 0.0151 | 0.0036 | 0.0130 | 0.0038 | 0.0140

(a) RME-GAN

(b) RadioUnet

(c) MBI

(d) Kriging

(e) RBF

(f) AE

Figure 7: Visualization Results of Different Methods

####  V-A3 Parameters of Learning Models

We utilize an 18-layer Unet for the generator with parameters shown in Table.
I, and a 4-layer discriminator shown as Fig. 3(b). In the first phase, the
weights of each loss term are set as $\lambda_{2}=10$, $\lambda_{3}=0.01$,
$\lambda_{4}=1$. In the second phase, we set the weights as $\lambda_{1}=10$,
$\lambda_{2}=1$, $\lambda_{3}=0.001$, $\lambda_{4}=0.0001$,
$\lambda_{5}=0.0001$ and $\lambda_{6}=0.0001$. Note that, although both phases
contains the $L_{\rm MSE}$, $\lambda_{2}$ is set to be smaller in the second
phase since we favor the downsampled radio map for detail correction.

###  V-B Overall Accuracy

We first present the results in different setups. To evaluate the performance
of the proposed RME-GAN, we compare with the MBI interpolation [4], RBF
interpolation [6], radioUnet [8], auto-encoder (AE) [10], deep AE [10],
traditional Unet, Kriging interpolation [17], and traditional cGAN [33]. The
Normalized Mean Square Error (NMSE) and Root Mean Square Error (RMSE) of the
reconstructed radio map are present in Table II.

In Setup 1, samples are uniformly distributed. The learning-based approaches
have better performance than the model-based and interpolation-based methods
since the environment of these regions are more complex and the radio map is
sensitive to shadowing effect. Our RME-GAN outperforms all the methods in
Setup 1. Especially in comparison with traditional cGAN, our two-phase
training process provides more information and display superior performance.
We also compare different methods over different ratio of observations in
Setup 1, shown in Fig. 6. Our proposed RME-GAN has the best performance
results, as Table II further illustrates.

In Setup 2 and Setup 3, observations are biased among different regions and
non-uniformly distributed, respectively. RME-GAN still shows the best results
owing to the design of different sampling strategies. The experimental results
demonstrate the power of the proposed resampling methods and RME-GAN in radio
map reconstruction.

Beyond the numerical results obtained, we also present the visualization
results of a radio map in Fig. 5(a). The reconstructed results from different
methods are shown in Fig. 7. Since interpolation-based methods, such as MBI
and RBF cannot capture the impact of environment efficiently, it fails to
reconstruct the radio map in both three setups. Kriging-based methods provide
a blurred construction of radio map as it estimates PSD from neighborhood
data. Learning-based methods generally display better reconstruction. Compared
to RadioUnet and AE, our RME-GAN constructs more accurate radio maps.
Especially in view of Setup 3 where the sparse observations are non-uniformly
distributed over the region, the proposed RME-GAN outperforms other existing
learning methods. This happens because our resampling strategy is able to
balance the samples over different locations. We also display the histogram of
the reconstructed radio map in Setup 1 from RME-GAN as shown in Fig. 8, where
RME-GAN is able to capture the distribution of PSD.

(a) Histogram in a training sample

(b) Histogram in a test sample

Figure 8: One Example of Histograms of PSD in Reconstructed Radio Maps Figure
9: Validation Performance over Different Phase

###  V-C Robustness

####  V-C1 Two-Phase Training

We first confirm the robustness of the two-phase training strategies. The
validation error over time is shown in Fig. 9. From the results, the error
drops when the training moves from phase 1 to phase 2, which indicates that
the second phase could better correct detailed errors for the propagation
template model from the first phase by leveraging the shadowing information.

####  V-C2 Zero-shot Learning

Zero-shot learning concerns the recognition of new concepts or new scenarios
when a learning model is trained with generic datasets [45]. To evaluate the
robustness and generalization of the proposed RME-GAN, we train the generator
model by using data from Setup 1, but we test the performance in the data from
Setup 2 and 3, respectively. The results are shown in Table III. From the
results, we see that our RME-GAN can still obtain a superior performance in
the test dataset with different distributions while other learning-based
methods all observe a significant performance drop.

TABLE III: Performance of Zero-shot Learning Model | Setup 2 | Setup 3
---|---|---
NMSE | RMSE | NMSE | RMSE
cGAN | 0.1613 | 0.0918 | 0.1984 | 0.1018
AE | 0.2051 | 0.1038 | 0.2259 | 0.1094
Deep AE | 0.1870 | 0.0991 | 0.1981 | 0.1019
Unet | 0.0082 | 0.0208 | 0.0254 | 0.0367
RadioUnet | 0.0052 | 0.0165 | 0.0148 | 0.0280
RME-GAN | 0.0042 | 0.0150 | 0.0067 | 0.0190
TABLE IV: Performance in Selected Areas Model | Red Area | Green Area | Blue Area
---|---|---|---
RMSE | NMSE | RMSE | NMSE | RMSE | NMSE
RBF | 0.0574 | 0.2794 | 0.0929 | 0.3317 | 0.0868 | 0.2803
Kriging | 0.0675 | 0.3867 | 0.0890 | 0.3041 | 0.0819 | 0.2491
AE | 0.0545 | 0.2518 | 0.0831 | 0.2651 | 0.0711 | 0.1879
MBI | 0.1087 | 0.9520 | 0.1368 | 0.7183 | 0.0952 | 0.3372
cGAN | 0.0449 | 0.1714 | 0.0723 | 0.2009 | 0.0466 | 0.0807
RadioUnet | 0.0164 | 0.0228 | 0.0250 | 0.0240 | 0.0210 | 0.0163
RME-GAN | 0.0151 | 0.0190 | 0.0221 | 0.0189 | 0.0172 | 0.0110
TABLE V: Outage Evaluation (Prediction Error) Model | Threshod = 5 | Threshod = 25
---|---|---
Samples = 1 | Average | Samples = 1 | Average
RBF | 0.1165 | 0.1265 | 0.2303 | 0.2349
MBI | 0.1169 | 0.1391 | 0.1599 | 0.2359
Kriging | 0.1870 | 0.0991 | 0.1381 | 0.1520
AE | 0.1168 | 0.1525 | 0.1498 | 0.1586
RadioUnet | 0.0085 | 0.0122 | 0.0144 | 0.0158
RME-GAN | 0.0074 | 0.0105 | 0.0120 | 0.0141

####  V-C3 Performance in Areas with Lower Energy

Compared to the areas with larger energy, the areas with lower power spectrum
are more sensitive to the environment. To validate the performance of RME-GAN
in areas with lower power, we select three areas shown as Fig. 5(b) and test
different methods in these three regions. The performance is shown in Table
IV, where RME-GAN outperforms other methods in all areas and display more
robustness in the low-power areas.

###  V-D Outage Fault Diagnosis

We perform outage analysis as a fault diagnosis. In some realistic
applications, e.g. outage detection, we usually only need a coarse resolution
of radio map instead of the fine-resolution one. To evaluate the power of the
proposed RME-GAN in reconstructing radio maps for such applications, we also
compare the accuracy of the reconstructed outage maps. In our definition, a
location experiences power outage if its PSD falls below a preset threshold.
Then, we can segment the radio map into an outage map with normal and outage
areas, i.e., a binary segmentation. We compare the outage maps for all methods
against two exemplary thresholds, i.e., 5 and 25. We measure performance by
using segmentation errors in Table V. Moreover, we present the results of one
randomly selected regions, and the average of 10 randomly selected regions.
Our proposed RME-GAN provides a more accurate outage map reconstruction
compared to other approaches, which demonstrates the practicability of RME-GAN
in network fault detection applications.

##  VI Conclusions

In this work, we propose a two-phase learning frameworks (RME-GAN) for radio
map estimation from sparse RF measurements. Our RME-GAN first estimates a
template to capture the global radio propagation patterns based model-based
interpolation. In the second phase it learns the shadowing effects via
geometric and frequency downsampling. The experimental results demonstrate the
performance advantages and the robustness of the proposed RME-GAN, together
with the two-phase training strategies in radio map estimation.

The rapid growth and expansion of IoT and 5G systems make efficient estimation
and utilization of radio map become increasingly important. Radiomap
estimation are important tools for resource allocation and network planning.
One promising future direction may be fault diagnosis based on the sparse RF
observations from deployed sensors or reports from a few users. Another
fruitful direction would be the integration of model-based approaches with
learning machines. Our future works plan to consider these and other related
directions.

## References

  * [1] Y. Deng et al., “Radio environment map construction using super-resolution imaging for intelligent transportation systems,” in IEEE Access, vol. 8, pp. 47272-47281, Mar. 2020.
  * [2] S. Bi, J. Lyu, Z. Ding and R. Zhang, “Engineering radio maps for wireless resource management,” in IEEE Wireless Communications, vol. 26, no. 2, pp. 133-141, Apr. 2019.
  * [3] S. Debroy, S. Bhattacharjee and M. Chatterjee, “Spectrum map and its application in resource management in cognitive radio networks,” in IEEE Transactions on Cognitive Communications and Networking, vol. 1, no. 4, pp. 406-419, Dec. 2015.
  * [4] M. Lee and D. Han, “Voronoi tessellation based interpolation method for wi-fi radio map construction,” in IEEE Communications Letters, vol. 16, no. 3, pp. 404-407, Mar. 2012.
  * [5] J. A. Bazerque, G. Mateos and G. B. Giannakis, “Group-lasso on splines for spectrum cartography,” in IEEE Transactions on Signal Processing, vol. 59, no. 10, pp. 4648-4663, Oct. 2011.
  * [6] J. Krumm and J. Platt, “Minimizing calibration effort for an indoor 802.11 device location measurement system,”, Microsoft Research, Nov. 2003.
  * [7] S. Kuo and Y. Tseng, “Discriminant minimization search for largescale rf-based localization systems,” in IEEE Transactions on Mobile Computing, vol. 10, no. 2, pp. 291-304, Feb. 2011.
  * [8] R. Levie, Ç. Yapar, G. Kutyniok and G. Caire, “Radiounet: fast radio map estimation with convolutional neural networks,” in IEEE Transactions on Wireless Communications, vol. 20, no. 6, pp. 4001-4015, Jun. 2021
  * [9] O. Ronneberger, P. Fischer, and T. Brox, “U-net: convolutional networks for biomedical image segmentation,” in International Conference on Medical image computing and computer-assisted intervention, Munich, Germany, Nov. 2015, pp. 234-241.
  * [10] Y. Teganya and D. Romero, “Deep completion autoencoders for radio map estimation,” in IEEE Transactions on Wireless Communications, vol. 21, no. 3, pp. 1710-1724, Mar. 2022.
  * [11] S. C. Beck, J. L. Turner, and O. Kovo, “High-resolution radio maps of wolf-rayet galaxies: optically thick h ii regions?.” The Astronomical Journal, vol. 120, no. 1, pp. 244, 2000.
  * [12] O. Esrafilian and D. Gesbert, “3d city map reconstruction from uav-based radio measurements,” GLOBECOM 2017, Singapore, Dec. 2017, pp. 1-6.
  * [13] J. Johansson, W. A. Hapsari, S. Kelley and G. Bodog, “Minimization of drive tests in 3GPP release 11,” in IEEE Communications Magazine, vol. 50, no. 11, pp. 36-43, Nov. 2012.
  * [14] M. Mirza, and S. Osindero,“Conditional generative adversarial nets,” arXiv preprint arXiv:1411.1784, 2014.
  * [15] X. Fu, N. D. Sidiropoulos, J. H. Tranter and W. K. Ma, “A factor analysis framework for power spectra separation and multiple emitter localization,” in IEEE Transactions on Signal Processing, vol. 63, no. 24, pp. 6581-6594, Dec.15, 2015.
  * [16] H. Braham, S. B. Jemaa, G. Fort, E. Moulines and B. Sayrac, “Fixed rank kriging for cellular coverage analysis,” in IEEE Transactions on Vehicular Technology, vol. 66, no. 5, pp. 4212-4222, May 2017.
  * [17] G. Boccolini, G. Hernández-Peñaloza and B. Beferull-Lozano, “Wireless sensor network for Spectrum Cartography based on Kriging interpolation,” 2012 IEEE 23rd International Symposium on Personal, Indoor and Mobile Radio Communications - (PIMRC), Sydney, NSW, Australia, Sep. 2012, pp. 1565-1570.
  * [18] A. E. C. Redondi, “Radio map interpolation using graph signal processing,” in IEEE Comm. Letters, vol. 22, no. 1, pp. 153-156, Jan. 2018.
  * [19] X. Du, X. Liao, M. Liu and Z. Gao, “CRCLoc: a crowdsourcing-based radio map construction method for wifi fingerprinting localization,” in IEEE Internet of Things Journal, vol. 9, no. 14, pp. 12364-12377, Jul. 2022.
  * [20] S. Zhang, T. Yu, J. Tivald, B. Choi, F. Ouyang, and Z. Ding, “Exemplar-based radio map reconstruction of missing areas using propagation priority,” arXiv preprint arXiv:2209.04566, 2022.
  * [21] M. S. Riaz, H. N. Qureshi, U. Masood, A. Rizwan, A. Abu-Dayya and A. Imran, “Deep learning-based framework for multi-fault diagnosis in self-healing cellular networks,” 2022 IEEE Wireless Communications and Networking Conference (WCNC), Austin, TX, USA, Apr. 2022, pp. 746-751.
  * [22] D. Schäufele, R. L. G. Cavalcante and S. Stanczak, “Tensor completion for radio map reconstruction using low rank and smoothness,” 2019 IEEE 20th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC), Cannes, France, Jul. 2019, pp. 1-5.
  * [23] C. Parera, Q. Liao, I. Malanchini, C. Tatino, A. E. C. Redondi and M. Cesana, “Transfer learning for tilt-dependent radio map prediction,” in IEEE Transactions on Cognitive Communications and Networking, vol. 6, no. 2, pp. 829-843, Jun. 2020.
  * [24] S. K. Vankayala, S. Kumar, I. Roy, D. Thirumulanathan, S. Yoon and I. S. Kanakaraj, “Radio map estimation using a generative adversarial network and related business aspects,” 2021 24th International Symposium on Wireless Personal Multimedia Communications (WPMC), Okayama, Japan, Feb. 2021, pp. 1-6.
  * [25] Chun-Hao Liu, Hun Chang, Taesuh Park; Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) Workshops, 2020, pp. 498-499.
  * [26] H. Zou et al., “Adversarial learning-enabled automatic wifi indoor radio map construction and adaptation with mobile robot,” in IEEE Internet of Things Journal, vol. 7, no. 8, pp. 6946-6954, Aug. 2020.
  * [27] X. Wang, X. Wang, S. Mao, J. Zhang, S. C. G. Periaswamy and J. Patton, “Indoor radio map construction and localization with deep gaussian processes,” in IEEE Internet of Things Journal, vol. 7, no. 11, pp. 11238-11249, Nov. 2020.
  * [28] Y. Zeng, X. Xu, S. Jin and R. Zhang, “Simultaneous navigation and radio mapping for cellular-connected uav with deep reinforcement learning,” in IEEE Transactions on Wireless Communications, vol. 20, no. 7, pp. 4205-4220, Jul. 2021.
  * [29] K. Suto, S. Bannai, K. Sato, K. Inage, K. Adachi and T. Fujii, “Image-driven spatial interpolation With deep learning for radio map construction,” in IEEE Wireless Communications Letters, vol. 10, no. 6.
  * [30] H. J. Bae and L. Choi, ”Large-scale indoor positioning using geomagnetic field with deep neural networks,” 2019 IEEE International Conference on Communications (ICC), Shanghai, China, May 2019, pp. 1-6.
  * [31] T. Imai, K. Kitao and M. Inomata, “Radio propagation prediction model using Cconvolutional neural networks by deep learning,” 2019 13th European Conference on Antennas and Propagation (EuCAP), Krakow, Poland, Jun. 2019, pp. 1-5.
  * [32] D. Romero and S. J. Kim, “Radio map estimation: a data-driven approach to spectrum cartography,” in IEEE Signal Processing Magazine, vol. 39, no. 6, pp. 53-72, Nov. 2022.
  * [33] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio. “Generative adversarial networks,” Communications of the ACM, vol. 63, pp. 139-144, Nov. 2020.
  * [34] D. Pathak, P. Krahenbuhl, J. Donahue, T. Darrell, and A. A. Efros, “Context encoders: Feature learning by inpainting,” Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, Las Vegas, Nevada, USA, 2016, pp. 2536-2544.
  * [35] Y. Liu, Z. Qin, T. Wan, and Z. Luo, “Auto-painter: cartoon image generation from sketch by using conditional Wasserstein generative adversarial networks,” Neurocomputing, vol. 311, pp. 78-87, Oct. 2018.
  * [36] P. Isola, J. Y. Zhu, T. Zhou, and A. A. Efros, “Image-to-image translation with conditional adversarial networks, in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Honolulu, Hawaii, USA, 2017, pp. 1125-1134.
  * [37] T. C. Wang, M. Y. Liu, J. Y. Zhu, A. Tao, J. Kautz, and B. Catanzaro, “High-resolution image synthesis and semantic manipulation with conditional gans,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, Salt Lake City, Utah, USA, 2018, pp. 8798-8807.
  * [38] S. Zhang, Q. Deng, and Z. Ding,“Multilayer graph spectral analysis for hyperspectral images,” EURASIP Journal on Advances in Signal Processing, vol. 1, no. 92 pp. 1-25, 2022.
  * [39] S. Zhang, Z. Ding and S. Cui, “Introducing hypergraph signal processing: theoretical foundation and practical applications,” in IEEE Internet of Things Journal, vol. 7, no. 1, pp. 639-660, Jan. 2020.
  * [40] H. V. Nguyen, and L. Bai, “Cosine similarity metric learning for face verification,” in Asian Conference on Computer Vision, Berlin, Heidelberg, 2010, pp. 709-720.
  * [41] Z. Wang, E. P. Simoncelli and A. C. Bovik, “Multiscale structural similarity for image quality assessment,” the Thrity-Seventh Asilomar Conference on Signals, Systems and Computers, Pacific Grove, CA, USA, 2003, pp. 1398-1402.
  * [42] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error measurement to structural similarity,” IEEE Trans. Image Processing, vol. 13, Jan. 2004.
  * [43] OpenStreetMap contributors. (2017). Planet Dump. [Online]. Available:https://planet.osm.org and https://www.openstreetmap.org.
  * [44] R. Hoppe, G. Wolfle, and U. Jakobus, “Wave propagation and radio network planning software WinProp added to the electromagnetic solver package FEKO,” in Proc. Int. Appl. Comput. Electromagn. Soc. Symp. Italy (ACES), Florence, Italy, Mar. 2017, pp. 1–2.
  * [45] B. R. Paredes, and P. Torr, “An embarrassingly simple approach to zero-shot learning,” in International Conference on Machine Learning, Lille, France, Jun. 2015, pp. 2152-2161.