# Feasibility study on solving the
Helmholtz equation in 3D with PINNs

Stefan Schoder
Institute of Fundamentals and Theory in Electical Engineering (IGTE)
Graz University of Technology
8010 Graz, Austria
<EMAIL_ADDRESS>
& Florian Kraxberger
Institute of Fundamentals and Theory in Electical Engineering (IGTE)
Graz University of Technology
8010 Graz, Austria
<EMAIL_ADDRESS>

###### Abstract

Room acoustic simulations at low frequencies often face significant
uncertainties of material parameters and boundary conditions due to absorbing
material. We discuss the application of Physics-Informed Neural Networks
(PINNs) to solve the (forward) Helmholtz equation in three dimensions (3D),
employing mini-batch stochastic gradient descent with periodic resampling
every 100 iterations for memory-efficient training. Addressing the
computational challenges posed by the extension of PINNs from 2D to 3D for
acoustics, DeepXDE is used for implementing the forward PINN. The proposed
numerical method is benchmarked against an analytical solution of a standing
wave field in 3D. The PINN results are also compared to the Finite Element
Method (FEM) solutions for a 3D wave field computed with openCFS. The
alignment between PINN-generated solutions and analytical/FEM solutions shows
the feasibility of PINNs modeling 3D acoustic applications for future inverse
problems, and validating the accuracy and reliability of the proposed
approach. Compared to FEM, establishing the PINN model took few hours (similar
to the setup of a FEM simulation), the training took 38h to 42.8h (which is
longer than the solution of the FEM simulation, which took 17min-19min), and
the inference took 0.05 seconds being more than 20,000 times faster than the
FEM benchmark openCFS using the same number of degrees of freedomwhen
producing the results. Thereby, the insight is gained that 3D acoustic wave
simulations in the frequency domain are feasible for forward PINNs and can
predict complex wave behaviors in real-world applications.

_K_ eywords PINNs $\cdot$ FEM $\cdot$ Helmholtz equation $\cdot$ DeepXDE
$\cdot$ openCFS $\cdot$ Acoustics $\cdot$ Waves $\cdot$ 3D $\cdot$ PDE $\cdot$
Room Acoustics $\cdot$ Time-Harmonic Wave Field $\cdot$ Absorber

##  1 Introduction

Room acoustics simulations aim towards improving the acoustic performance of
interior spaces by meeting specific criteria and predicting acoustic quality
in advance of the construction phase. However, numerous challenges have to be
overcome, including uncertain boundary conditions and the simulation of low-
frequent wave phenomena [1, 2]. A comprehensive examination underscores the
current obstacles, with a particular emphasis on the persisting challenge of
errors arising from uncertain input parameters. A survey identified that about
two-thirds of the acousticians judge uncertainties in input parameters as the
foremost obstacle hindering more accurate simulation results [3]. In [4],
simulations with the use of absorption coefficients from textbook tables
showed inadequate prediction capabilities of reverberation time. These
insights underscore the critical importance of using correct material data for
the accurate modeling of room acoustics simulations, manifesting the interest
in using forward physics-informed neural networks (PINNs) and inverse PINNs in
three-dimensional (3D) space [1]. The computation of the time-harmonic
acoustic wave field in large 3D geometries represents a challenging frontier
for numerical methods. Especially, the solutions to the Helmholtz equation in
3D are essential for understanding low-frequent room acoustic scenarios.
Traditional numerical methods, such as Finite Element Methods (FEM), have
provided valuable insights [5, 6], but their computational approach is
limited, and the inverse estimation of parameters is challenging. In this
context, PINNs show a promising supplement, leveraging the power of
optimization and machine learning to tackle both the forward propagation of
acoustics and the inverse estimation of material parameters of 3D acoustic
fields.

This paper addresses the challenge associated by extending PINNs to solve the
Helmholtz equation in 3D, presenting a methodological breakthrough facilitated
by the integration of state-of-the-art mini-batch stochastic gradient descent
optimization with periodic shuffling of the training data points. The PINN
investigation is conducted within the DeepXDE framework [7] and compared to a
continuously tested and validated [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 2]
FEM solver openCFS [18]. Physics-Informed Neural Networks offer a unique
advantage, using the same optimization routines when establishing the forward
PINN and the inverse PINN. By incorporating domain-specific knowledge and
physical principles into the neural network architecture, PINNs have
demonstrated a capacity to learn and predict complex solutions to partial
differential equations [19, 20]. This study builds upon the success of PINNs
in 1D and 2D acoustic simulations [1, 7, 21, 22, 23, 24] (see Tab. 1,
extending their application to the more challenging 3D acoustics problems. The
goal is to model and predict wave phenomena in 3D space accurately.

Table 1: Literature overview on noteworthy public development projects solving acoustic fields. Type |  Name |  URL |  Ref.
---|---|---|---
Article, GIT |  DeepXDE: A Deep Learning Library for Solving Differential Equations |  GIT |  [7]
Article, GIT |  Solving the frequency-domain acoustic VTI wave equation using physics-informed neural networks |  GIT |  [21]
Article |  Hyper-parameter tuning of physics-informed neural networks: Application to Helmholtz problems |  |  [22]
Preprint |  FO-PINNs: A First-Order formulation for Physics Informed Neural Networks |  |  [23]
Article |  Helmholtz-equation solution in nonsmooth media by a physics-informed neural network incorporating quadratic terms and a perfectly matching layer condition |  |  [24]
Preprint |  PINNs-TF2: Fast and User-Friendly Physics-Informed Neural Networks in TensorFlow V2 |  GIT |  [25]
Preprint |  Physics-Informed Neural Networks for Acoustic Boundary Admittance Estimation |  Preprint |  [26]

The article is structured as follows: Sec. 2 introduces the theoretical
background of PINNs as self-supervised learning. Section 3 presents the
numerical examples of 3D room-like geometry with different boundary
conditions, differential equation forcing, and material parameters. This
section also shows the training results and the comparison to the analytic
solution and FEM solution by error analysis. Each numerical example is
accompanied by a discussion of the results obtained. Finally, Sec. 4
summarizes the key findings and provides a roadmap for future research
directions.

##  2 Physics-informed Neural Networks for the Helmholtz Equation in Acoustics

PINNs represent a cutting-edge fusion of neural networks and equation-based
principles of physics, offering a bundle of possibilities for solving complex
partial differential equations (PDEs) and simulating physical systems [19].
This innovative approach delegates implementation-associated challenges
encountered in traditional computational methods, such as Finite Element
Methods (FEM) or Finite Difference Methods, by seamlessly integrating
computational domain-specific and mathematical equation-specific knowledge
directly into the neural network architecture. Specifically, Figure 1
schematically shows how such a feed-forward (multi-layer perceptron [27, 28])
fully connected layer network is tailored to model the inhomogeneous Helmholtz
equation

$(\Delta+k^{2})p(\omega,\bm{x})=f(\omega,\bm{x})$ (1)

where $p\in\mathbb{C}^{3}$ is the time-harmonic acoustic pressure,
${k}=\omega/c$ is the wavenumber, $\omega$ the angular frequency, $c$ the
speed of sound, $f$ the forcing term, and $\Delta$ is the Laplacian operator.

To predict $\hat{p}(\omega,\bm{x})$ with a PINN, let
$\mathcal{N}^{L}(\mathbf{x}):\mathbb{R}^{d_{\text{in
}}}\rightarrow\mathbb{R}^{d_{\text{out }}}$ be a $L$-feed-forward neural
network (FNN) (with $(L-1)$-hidden neural network layers), with $N_{\ell}$
neurons in the $\ell$-th layer $\left(N_{0}=d_{\text{in }},N_{L}=d_{\text{out
}}\right)$. The output at layer of a feed-forward neural network is calculated
by repeated nonlinear activation-function $\sigma$ weighted tensor products of
(i) inputs and (ii) bias vector and weight matrix (e.g. in the $\ell$-th layer
by $\bm{W}^{\ell}\in\mathbb{R}^{N_{\ell}\times N_{\ell-1}}$ and
$\mathbf{b}^{\ell}\in\mathbb{R}^{N_{\ell}}$, respectively), such that the FNN
111Note that all the network parameters depend on the solution frequency
$\omega$. When the frequency of the problem is changing, the FNN must be
retrained. This is similar to the FEM, where the solver must construct and
invert the system matrices for different frequencies separately. is
recursively defined by

input layer:
$\displaystyle\mathcal{N}^{0}(\mathbf{x})=\mathbf{x}\in\mathbb{R}^{d_{\mathrm{in}}}$
(2) hidden layers:
$\displaystyle\mathcal{N}^{\ell}(\mathbf{x})=\sigma\left(\bm{W}^{\ell}\mathcal{N}^{\ell-1}(\mathbf{x})+\bm{b}^{\ell}\right)\in\mathbb{R}^{N_{\ell}},\quad\text{
for }\quad 1\leq\ell\leq L-1,$ output layer:
$\displaystyle\mathcal{N}^{L}(\mathbf{x})=\bm{W}^{L}\mathcal{N}^{L-1}(\mathbf{x})+\bm{b}^{L}\in\mathbb{R}^{d_{\text{out
}}};$

In the case of PINNs in 3D, $\bm{x}=(x,y,z)$ at the input layer, where
$(x,y,z)$ denotes the space directions of one point. The output prediction at
the output layer $L$ is calculated by

$\displaystyle\hat{p}=$
$\displaystyle\bm{W}^{L}\mathcal{N}^{L-1}(\mathbf{x})+\bm{b}^{L}=\bm{W}^{L}\sigma\left(\bm{W}^{L-1}\mathcal{N}^{L-2}(\mathbf{x})+\bm{b}^{L-1}\right)+\bm{b}^{L}=\dots$
(3) $\displaystyle=$
$\displaystyle\mathcal{N}^{L}\circ\mathcal{N}^{L-1}\circ\dots\circ\mathcal{N}^{0}\circ\bm{x}=\mathcal{N}(\bm{\theta},\bm{x})\,,$
(4)

with $\bm{g}$ being a series of mapping functions. The network transformative
operations can be summarized by a nonlinear operation
$\mathcal{N}(\bm{\theta},\bm{x})$ and a structure $\bm{\theta}$ collecting all
weights and biases $\bm{\theta}=\\{\bm{W}^{l},\bm{b}^{l}\\}_{1\leq l\leq L}$.

Figure 1: PINNs applied to Helmholtz equation.

During supervised learning, the neural network weights will be learned by
labeled training data $\mathcal{D}=\\{\bm{x}_{i},p_{i}\\}_{i=1}^{N_{L}}$ of
$N_{L}$ samples, and unlabeled training data
$\mathcal{T}=\\{\bm{x}_{i}\\}_{i=1}^{N_{U}}$ of $N_{U}$ samples, to
approximate a continuous function mapping $g(\bm{x})=p$. $\mathcal{T}$
comprises two sets, the points in the domain $\mathcal{T}_{f}$ of $N_{PDE}$
samples and the points on the boundary $\mathcal{T}_{b}$ of $N_{BC}$ samples.
The neural network operator $\hat{p}=\mathcal{N}(\bm{\theta},\bm{x})$ can be
interpreted as an approximation of $g(\bm{x})=p$, regarding optimality of the
discrepancy (error) described by the cost functional (loss function), e.g.
based on mean squared error (MSE) for the data loss

$\mathcal{L}_{\mathrm{data}}(\bm{\theta},\mathcal{D})=\frac{1}{N_{L}}\sum_{i=1}^{N}\parallel
p_{i}-\hat{p}_{i}\parallel^{2}$ (5)

with $\|\cdot\|$ being the $L_{2}$-norm. In general, the data loss can be
combined with other types of loss (e.g., boundary loss or PDE loss) as
described in sec. 3, resulting in a total loss $\mathcal{L}_{\mathrm{total}}$.
According to the minimization problem

$\hat{\bm{\theta}}=\underset{\bm{\theta}}{\arg\min}\,(\lambda_{\mathrm{data}}\mathcal{L}_{\mathrm{data}}(\bm{\theta},\mathcal{D})+\lambda_{\mathrm{pde}}\mathcal{L}_{\mathrm{pde}}(\bm{\theta},\mathcal{T}),$
(6)

the weights and biases $\bm{\theta}$ are adjusted to minimize the loss
function $\mathcal{L}$. The converged set of weights and biases $\bm{\theta}$
is the optimal prediction of these parameters, denoted by $\hat{\bm{\theta}}$.
Minimization can be achieved by gradient-based optimization algorithms, like
Adam [29] or LBFGS algorithm [30]. Using the gradient of the loss function
$\nabla_{\bm{\theta}}\mathcal{L}$, the weights and biases are adapted
iteratively to decrease $\mathcal{L}$, typically using backpropagation [31]
and automatic differentiation [32]. A standard gradient descent algorithm
updates the parameters such that

$\bm{\theta}_{i+1}=\bm{\theta}_{i}-\alpha_{i}\nabla_{\bm{\theta}_{i}}\mathcal{L}\,,$
(7)

with $\alpha$ representing the learning rate and $i$ the epoch number.

The Helmholtz equation (1) is a second-order linear partial differential
equation of hyperbolic type governing the behavior of acoustic waves in the
frequency domain. It is derived from the linear scalar wave equation in the
time domain through separation of variables and an assumption of time-harmonic
dependence $\tilde{p}(\bm{x},t)=p(\bm{x})e^{-\mathrm{i}\omega t}$. In addition
to (1), for the inhomogeneous Helmholtz equation on a restricted domain
$\Omega$, forcing and boundary conditions (e.g., Dirichlet boundary condition
at $\Gamma_{D}$ and Neumann boundary $\Gamma_{N}$, where the two surfaces form
the domain boundary denoted as $\partial\Omega=\Gamma_{D}\cup\Gamma_{N}$) have
to be considered, forming a well-posed system of equations

$\displaystyle\Delta p(\bm{x})+k^{2}p(\bm{x})$
$\displaystyle=f(x),\quad\bm{x}\in\Omega\setminus\partial\Omega\,,$ (8a)
$\displaystyle p(\bm{x})$ $\displaystyle=0,\quad\bm{x}\in\Gamma_{D}\,,$ (8b)
$\displaystyle\nabla p(\bm{x})\cdot\bm{n}$
$\displaystyle=0,\quad\bm{x}\in\Gamma_{N}\,.$ (8c)

Here, $f(\bm{x})$ represents the known volume forcing term of the
inhomogeneous Helmholtz equation and $\bm{n}$ is the outward pointing normal
vector of $\Gamma_{N}$. The solution to the Helmholtz equation yields the
acoustic pressure $p(\bm{x})$, typically expressed as a complex number, at a
specific location $\bm{x}$. To incorporate the inherent physical knowledge
into the neural network, the residual of the Helmholtz equation
$r_{\mathrm{PDE}}$ is integrated into the loss function. The residual of the
PDE is

$\begin{split}r_{\mathrm{PDE}}(\bm{x})&=\Delta\hat{p}(\bm{x})+k^{2}\hat{p}(\bm{x})-f(\bm{x})\qquad\forall\bm{x}\in
D_{\mathrm{PDE}}\\\
\mathcal{L}_{\mathrm{PDE}}(\bm{\theta},\mathcal{T}_{f})&=\frac{1}{N_{\mathrm{PDE}}}\sum_{i=1}^{N_{\mathrm{PDE}}}\parallel
r_{\mathrm{PDE}}(\bm{x})\parallel^{2}\end{split}$ (9)

and it is evaluated at randomly sampled collocation points $\bm{x}$ within the
computational domain $\Omega$, forming the dataset
$\mathcal{T}_{f}=\\{\bm{x}_{i}\\}_{i=1}^{N_{\mathrm{PDE}}}$, where
$N_{\mathrm{PDE}}$ is the number of sample points for the PDE residual. In
addition, the residuals of the boundary conditions have to be integrated into
the loss function too. The Dirichlet boundary condition residual
$r_{\mathrm{DBC}}$

$\begin{split}r_{\mathrm{DBC}}&=p(\bm{x})\qquad\forall\bm{x}\in
D_{\mathrm{DBC}}\,,\\\
\mathcal{L}_{\mathrm{DBC}}(\bm{\theta},\mathcal{T}_{b1})&=\frac{1}{N_{\mathrm{DBC}}}\sum_{i=1}^{N_{\mathrm{DBC}}}\parallel
r_{\mathrm{DBC}}(\bm{x})\parallel^{2}\end{split}$ (10)

is sampled at the respective sampling points on the Dirichlet boundary
$\Gamma_{D}$ forming the set
$\mathcal{T}_{b1}=\\{\bm{x}_{i}\\}_{i=1}^{N_{\mathrm{DBC}}}$, where
$N_{\mathrm{DBC}}$ is the number of sample points for the Dirichlet boundary
$\Gamma_{D}$. Similarly, the Neumann boundary condition residual
$r_{\mathrm{NBC}}$ is

$\begin{split}r_{\mathrm{NBC}}&=\nabla
p(\bm{x})\cdot\bm{n}\qquad\forall\bm{x}\in D_{\mathrm{NBC}}\,\\\
\mathcal{L}_{\mathrm{NBC}}(\bm{\theta},\mathcal{T}_{b2})&=\frac{1}{N_{\mathrm{NBC}}}\sum_{i=1}^{N_{\mathrm{NBC}}}\parallel
r_{\mathrm{NBC}}(\bm{x})\parallel^{2}\,,\end{split}$ (11)

which is sampled at the respective sampling points on the Neumann boundary
$\Gamma_{N}$ forming the set
$\mathcal{T}_{b2}=\\{\bm{x}_{i}\\}_{i=1}^{N_{\mathrm{NBC}}}$, where
$N_{\mathrm{NBC}}$ is the number of sample points for the Neumann boundary
$\Gamma_{B}$. In the training process, incorporating the BCs involves adding
the residual of the BCs directly to the loss function as a mean squared error,
corresponding to the forward solution. This results in a total loss

$\mathcal{L}_{\mathrm{total}}=\lambda_{\mathrm{data}}\mathcal{L}_{\mathrm{data}}+\lambda_{\mathrm{PDE}}\mathcal{L}_{\mathrm{PDE}}+\lambda_{\mathrm{DBC}}\mathcal{L}_{\mathrm{DBC}}+\lambda_{\mathrm{NBC}}\mathcal{L}_{\mathrm{NBC}}\,,$
(12)

where $\lambda_{\mathrm{data}}$, $\lambda_{\mathrm{PDE}}$,
$\lambda_{\mathrm{DBC}}$, and $\lambda_{\mathrm{NBC}}$ are weighting factors
for the individual loss terms. If not defined otherwise
$\lambda_{\mathrm{PDE}}=1$ for all used FNN.

##  3 Application

In this section, basic room application examples in 2D and 3D are presented to
highlight the capabilities of PINNs when approximating solutions of the
Helmholtz equation.

###  3.1 Example 1a - Analytic solution in 2D, Dirichlet

The example demonstrates222Details of this initial example can be found here
https://deepxde.readthedocs.io/en/latest/demos/pinn_forward/helmholtz.2d.dirichlet.html.
the application of PINNs to solve the Helmholtz equation in a 2D domain with
Dirichlet boundary conditions. The governing equation is eq. (1) with forcing
of

$f(x,y)=k^{2}\sin(kx)\sin(ky)\,,$ (13)

with the wave number of $k=2\pi/\lambda$, wavelength $\lambda=1/2$, and the
solution is sought in a 2D spatial domain $\Omega=[0,1]^{2}$. This leads to
the analytic solution

$p(x,y)=\sin(kx)\sin(ky)\,.$ (14)

The training process involves formulating a comprehensive loss function that
incorporates the Helmholtz equation. For the Helmholtz equation, the loss term
is defined in eq. (9). A homogeneous Dirichlet boundary condition $p(x,y)=0$
with $(x,y)\in\partial\Omega$ is enforced by a transform of the neural network
solution $\mathcal{N}(x,y)$ to the acoustic pressure

$\tilde{p}(x,y)=x(x-1)y(y-1)\mathcal{N}(x,y)\,.$ (15)

The resulting total loss function for this example is equivalent to eq. (12)
with
$\lambda_{\mathrm{data}}=\lambda_{\mathrm{NBC}}=\lambda_{\mathrm{DBC}}=0$, and
therefore reduces to

$\mathcal{L}_{\mathrm{1a}}=\mathcal{L}_{\mathrm{PDE}}\,.$ (16)

The PINN is set up using DeepXDE (with the PyTorch backend), and the trained
network is capable of predicting the acoustic pressure $\tilde{p}$ within the
specified 2D domain. Ten random collocation points per wavelength along each
direction for training and 30 for testing are defined. A fully connected
neural network of four layers depth (three hidden layers) and layer width of
150 neurons is used. A sinus activation function is used with Glorot uniform
bias and weights initialization. The network is trained over 5000 iterations
by the ADAM optimizer with a learning rate of 0.001, resulting in a test loss
of 0.00571.444The defined hyperparameters can be optimized using Ray Tune or
scikit-
optimize333https://deepxde.readthedocs.io/en/latest/demos/pinn_forward/helmholtz.2d.dirichlet.hpo.html.

To enforce the boundary conditions via the loss function using
$\lambda_{\mathrm{DBC}}=100$, 15000 iterations have been made, resulting in a
test loss of 0.0761. Figure 2 shows the exact (analytical) solution, the PINN
solution with constraint DBC, and the PINN solution with DBC modeled as an
additional loss term. The figure shows the acoustic pressure of the whole
domain.

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex1-2D/plot_amplitude_exact.png} \put(-4.0,0.0){0}
\put(1.0,-6.0){0} \put(-4.0,84.0){1} \put(85.0,-6.0){1} \put(99.0,1.0){-1}
\put(99.0,42.0){~{}0} \put(99.0,84.0){~{}1} \end{overpic} (a) Exact solution.

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex1-2D/plot_amplitude.png} \put(1.0,-6.0){0}
\put(85.0,-6.0){1} \put(99.0,1.0){-1} \put(99.0,42.0){~{}0}
\put(99.0,84.0){~{}1} \end{overpic} (b) PINN, constraint DBC.

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex1-2D/plot_amplitude_weak.png} \put(1.0,-6.0){0}
\put(85.0,-6.0){1} \put(99.0,1.0){-1} \put(99.0,42.0){~{}0}
\put(99.0,84.0){~{}1} \end{overpic} (c) PINN, DBC loss.

Figure 2: The real part of the acoustic pressure.

###  3.2 Example 1b - Analytic solution in 3D, Dirichlet

Consistently with the previous example in 3D the computational domain is
denoted by $\Omega=[0,1]^{3}$, the application of PINNs to solve the Helmholtz
equation is extended using Neumann boundary conditions. The governing equation
is eq. (1) with forcing

$f(x,y,z)=2k^{2}\cos(kx)\cos(ky)\cos(kz)\,,$ (17)

with a wave number of $k=2\pi/\lambda$, wave length $\lambda=1/2$. This leads
to an analytic solution of

$p(x,y,z)=\cos(kx)\cos(ky)\cos(kz)\,.$ (18)

For the Helmholtz equation in 3D, the loss term is defined by eq. (9). A
homogeneous Dirichlet boundary condition $p(x,y,z)=0$ with
$(x,y,z)\in\partial\Omega$ is enforced by a transform of the neural network
solution $\mathcal{N}(x,y,z)$ to the acoustic pressure

$\tilde{p}(x,y,z)=x(x-1)y(y-1)z(z-1)\mathcal{N}(x,y,z)\,.$ (19)

The resulting total loss function for this example is equivalent to eq. (12)
with
$\lambda_{\mathrm{data}}=\lambda_{\mathrm{NBC}}=\lambda_{\mathrm{DBC}}=0$, and
therefore reduces to the one from the previous example (16). Similarly to the
example shown in sec. 3.1, the PINN is set up using DeepXDE (PyTorch backend)
with an updated input layer for the additional space domain.

The trained network is capable of predicting the acoustic pressure $\tilde{p}$
within the specified 3D domain, with the same collocation points resolution as
defined in the previous example. A fully connected neural network of four
layers depth (three hidden layers) and layer width of 250 is defined. Instead
of the gradient descent using all collocation points, a batch-gradient descent
algorithm is used for RAM efficient optimization of the neural network weights
with periodic random resampling of the training points after 100 iterations. A
sinus activation function is used with Glorot uniform bias and weights
initialization. The network is trained over 10000 iterations using the ADAM
optimizer with a learning rate of 0.001, resulting in a test loss of 0.142.

In the case of enforcing the boundary conditions via the loss function using
$\lambda_{\mathrm{DBC}}=100$, 20000 iterations (layer-width 180 have been
made, resulting in a test loss of 0.0152. Figure 3 shows the exact solution,
the PINN solution with constraint DBC, and the PINN solution with DBC modeled
as an additional loss term. The figure shows the acoustic pressure in a cut
through the whole domain at $z=0.125$. It is observed that the number of
iterations required rises with the number of loss terms and rises with the
spatial dimensions.

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex1-3D/plot_amplitude_exact.png} \put(-4.0,0.0){0}
\put(1.0,-6.0){0} \put(-4.0,84.0){1} \put(85.0,-6.0){1} \put(99.0,1.0){-1}
\put(99.0,42.0){~{}0} \put(99.0,84.0){~{}1} \end{overpic} (a) Exact solution.

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex1-3D/plot_amplitude.png} \put(1.0,-6.0){0}
\put(85.0,-6.0){1} \put(99.0,1.0){-1} \put(99.0,42.0){~{}0}
\put(99.0,84.0){~{}1} \end{overpic} (b) PINN, constraint DBC.

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex1-3D/plot_amplitude_85.png} \put(1.0,-6.0){0}
\put(85.0,-6.0){1} \put(99.0,1.0){-1} \put(99.0,42.0){~{}0}
\put(99.0,84.0){~{}1} \end{overpic} (c) PINN, DBC loss.

Figure 3: The real part of the acoustic pressure at $z=0.125$.

Until now, no feasible boundary conditions for room acoustic applications have
been used. In room acoustic simulations, it is justified to use homogeneous
Neumann boundary conditions as an approximation of the impedance jump between
air and the room wall for low frequencies [2]. As a next step, the loss-term
modeled homogeneous Dirichlet boundary condition (sound soft) will be replaced
by a homogeneous Neumann boundary condition (sound hard wall).

###  3.3 Example 2a - Analytic solution in 3D, Room Acoustics

Following the previous example, the application of PINNs to solve the
Helmholtz equation is extended to a 3D domain $\Omega=[0,1]^{3}$ with enforced
Dirichlet boundary conditions. The governing equation is eq. (1) with forcing
of

$f(x,y,z)=2k^{2}\sin(kx)\sin(ky)\sin(kz)\,,$ (20)

with a wave number of $k=2\pi/\lambda$ and wave length $\lambda=1/2$. This
leads to an analytic solution of

$p(x,y,z)=\sin(kx)\sin(ky)\sin(kz)\,.$ (21)

A homogeneous Neumann boundary condition $\nabla p(x,y,z)\cdot\bm{n}=0$ with
$(x,y,z)\in\partial\Omega$ is considered by the loss function. The loss term
is defined in eq. (12), with
$\lambda_{\mathrm{data}}=\lambda_{\mathrm{DBC}}=0$,
$\lambda_{\mathrm{NBC}}=5$, and $\lambda_{\mathrm{PDE}}=1$. The PINN is set up
using DeepXDE (PyTorch backend) with a layer-width of 180555The model
coincides with the loss-term modeled Dirichlet boundary condition Helmholtz
PINN from the previous example and can be seen as an evolution of this model..
The network is trained over 30000 iterations by the ADAM optimizer with a
learning rate of 0.001, resulting in a test loss of 0.125.

Figure 4 shows the exact solution, the PINN solution with constraint DBC, and
the PINN solution with DBC modeled as an additional loss term. The figure
shows the acoustic pressure in a cut through the whole domain at $z=0.125$. It
is observed that the number of iterations required rises with the use of a
Neumann boundary condition. Furthermore, iterating the final ADAM optimized
network parameters with an L-BFGS optimizer for 15000 iterations brings down
the test loss to 0.0142.

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex2a/plot_amplitude_exact.png} \put(-4.0,0.0){0}
\put(1.0,-6.0){0} \put(-4.0,84.0){1} \put(85.0,-6.0){1} \put(99.0,1.0){-1}
\put(99.0,42.0){~{}0} \put(99.0,84.0){~{}1} \end{overpic} (a) Exact solution.

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex2a/plot_amplitude.png} \put(1.0,-6.0){0}
\put(85.0,-6.0){1} \put(99.0,1.0){-1} \put(99.0,42.0){~{}0}
\put(99.0,84.0){~{}1} \end{overpic} (b) PINN ADAM.

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex2b/plot_amplitude_noSigma.png} \put(1.0,-6.0){0}
\put(99.0,1.0){-1} \put(99.0,42.0){~{}0} \put(99.0,84.0){~{}1} \end{overpic}
(c) PINN ADAM & L-BFGS.

Figure 4: The real part of the acoustic pressure at $z=0$ with NBC.

Another hyperparameter test was trying to predict both the real and the
imaginary part of the acoustic pressure $p$, wherein the imaginary part is
zero in the following example. As a consequence, when using the PINN for
obtaining a complex acoustic pressure, it is suggested to train a network for
the real and the imaginary part separately.

###  3.4 Example 2b - Comparison to FEM in 3D, Room Acoustics

As a follow-up to the previous 3D room acoustic example, the source is modeled
more realistically by a point source in the center of the room modeled by

$f(x,y,z)=2k^{2}\cos(kx)\cos(ky)\cos(kz)\mathrm{e}^{-\frac{(x-0.5)^{2}+(y-0.5)^{2}+(z-0.5)^{2}}{2\sigma^{2}}}\,.$
(22)

The results for this excitation, are studied for several parameters $\sigma$,
being a measure for the sharpness of the source distribution. For
$\sigma\rightarrow 0$, the source is very localized. In the current form, the
source decreases it strength with $\sigma\rightarrow 0$. Until a parameter
$\sigma>0.1$, the original resolution of the field quantities was well enough
to resolve the source distribution, leading to a condition $\lambda<6\sigma$.
As soon as this condition, does not hold, the resolution of the training data
is too rough for the PINN to be accurate. For smaller $\sigma<0.1$, additional
adaptive refinement was introduced, enriching the number of training points
(recursively by five, for 10
iterations666https://deepxde.readthedocs.io/en/latest/demos/pinn_forward/burgers.rar.html)
in areas where the residual $r_{\mathrm{PDE}}$ were largest. After the
enriching of training points, 3000 iterations were performed with the ADAM
optimization following by 150000 iterations using the L-BFGS. With this
strategy, $0.1>\sigma>0.02$ were predicted reasonable and a first step to
resolve spatial and amplitude multi-scales of the problem.

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex2b/plot_amplitude_noSigma.png} \put(-4.0,0.0){0}
\put(1.0,-6.0){0} \put(-4.0,84.0){1} \put(85.0,-6.0){1} \put(99.0,1.0){-1}
\put(99.0,42.0){~{}0} \put(99.0,84.0){~{}1} \end{overpic} (a)
$\sigma\rightarrow\infty$

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex2b/plot_amplitude_100.0.png} \put(1.0,-6.0){0}
\put(85.0,-6.0){1} \put(99.0,1.0){-1} \put(99.0,42.0){~{}0}
\put(99.0,84.0){~{}1} \end{overpic} (b) $\sigma=100$

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex2b/plot_amplitude_10.0.png} \put(1.0,-6.0){0}
\put(85.0,-6.0){1} \put(99.0,1.0){-1} \put(99.0,42.0){~{}0}
\put(99.0,84.0){~{}1} \end{overpic} (c) $\sigma=10$

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex2b/plot_amplitude_1.0.png} \put(-4.0,0.0){0}
\put(1.0,-6.0){0} \put(-4.0,84.0){1} \put(85.0,-6.0){1} \put(99.0,1.0){-1}
\put(99.0,42.0){~{}0} \put(99.0,84.0){~{}1} \end{overpic} (d) $\sigma=1$

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex2b/plot_amplitude_0.1.png} \put(1.0,-6.0){0}
\put(85.0,-6.0){1} \put(99.0,1.0){-1} \put(99.0,42.0){~{}0}
\put(99.0,84.0){~{}1} \end{overpic} (e) $\sigma=0.1$

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex2b/plot_amplitude_0.05.png} \put(1.0,-6.0){0}
\put(85.0,-6.0){1} \put(99.0,1.0){-1} \put(99.0,42.0){~{}0}
\put(99.0,84.0){~{}1} \end{overpic} (f) $\sigma=0.05$

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex2b/plot_amplitude_0.04.png} \put(-4.0,0.0){0}
\put(1.0,-6.0){0} \put(-4.0,84.0){1} \put(85.0,-6.0){1} \put(99.0,1.0){-1}
\put(99.0,42.0){~{}0} \put(99.0,84.0){~{}1} \end{overpic} (g) $\sigma=0.04$

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex2b/plot_amplitude_0.03.png} \put(1.0,-6.0){0}
\put(85.0,-6.0){1} \put(99.0,1.0){-1} \put(99.0,42.0){~{}0}
\put(99.0,84.0){~{}1} \end{overpic} (h) $\sigma=0.03$

\begin{overpic}[width=121.41306pt,trim=71.13188pt 31.2982pt 82.51282pt
0.0pt,clip]{media/Ex2b/plot_amplitude_0.02.png} \put(1.0,-6.0){0}
\put(85.0,-6.0){1} \put(99.0,1.0){-1} \put(99.0,42.0){~{}0}
\put(99.0,84.0){~{}1} \end{overpic} (i) $\sigma=0.02$

Figure 5: The real part of the acoustic pressure at $z=0.5$.

The FEM reference solution based on the same setup was carried out in openCFS
[18]. Thereby, a sparse and a fine mesh are used as depicted in fig. 6. The
sparse mesh discretizes the cube’s edge by 20 linear elements (approximately
10 degrees of freedom per wave length) , resulting in a total number of 8000
elements, and the fine mesh discretizes the cube’s edge by 80 linear elements
(approximately 40 degrees of freedom per wave length), resulting in 512000
elements of the computational domain.

(a) Sparse FEM mesh (FEM-sparse).

(b) Fine FEM mesh (FEM-fine).

Figure 6: Cross-section through the sparse and fine meshes used for the FEM
computations at $z=0.5$.

A relative error measure $Err_{\mathrm{rel}}$ is evaluated, such that

$Err_{\mathrm{rel}}=\sqrt{\frac{\sum_{i=1}^{N_{\mathrm{pts}}}\left(p_{i,\mathrm{ref}}-p_{i}\right)^{2}}{\sum_{i=1}^{N_{\mathrm{pts}}}\left(p_{i,\mathrm{ref}}\right)^{2}}}\,,$
(23)

where $N_{\mathrm{pts}}$ is the number of evaluation points,
$p_{i,\mathrm{ref}}$ is the reference pressure and $p_{i}$ is the pressure for
which the agreement against the reference pressure should be quantified. The
$N_{\mathrm{pts}}=10000$ evaluation points have been placed to the locations
where the PINN prediction is evaluated, i.e., a grid of
$$10\text{\,}\frac{\mathrm{cm}}{}$\times$10\text{\,}\frac{\mathrm{cm}}{}$$
across the computational domain at $z=$0.5\text{\,}\frac{\mathrm{m}}{}$$. To
interpolate the FEM solution to this grid, the FE basis functions are
evaluated at the according location in the respective element, as implemented
in openCFS [18]. The errors are denoted as follows:
$Err_{\mathrm{rel}}^{\text{PINN,FEMsparse}}$ is the error between the PINN
solution and the FEM solution obtained with the sparse mesh,
$Err_{\mathrm{rel}}^{\text{PINN,FEMfine}}$ is the error between the PINN
solution and the FEM solution obtained with the fine mesh, and
$Err_{\mathrm{rel}}^{\text{FEM}}$ is the error between the FEM solutions
obtained with fine and sparse meshes.

Table 2: Performance comparison between PINN and FEM. For FEM-sparse and FEM-fine the given durations are the CPU-hours, and for PINN the duration denotes the time on one GPU. $\sigma$ | Duration (mm:ss) | Error
---|---|---
PINN-training | PINN-prediction | FEM-sparse | FEM-fine | $Err_{\mathrm{rel}}^{\text{PINN,FEMsparse}}$ | $Err_{\mathrm{rel}}^{\text{PINN,FEMfine}}$ | $Err_{\mathrm{rel}}^{\text{FEM}}$
$0.1$ | 2374:40.4 | 00:00.05 | 00:16.8 | 17:44.8 | $0.5754$ | $0.9718$ | $0.9275$
$1.0$ | 2512:01.8 | 00:00.05 | 00:17.6 | 17:18.8 | $0.0454$ | $0.0997$ | $0.0924$
$10.0$ | 2280:20.0 | 00:00.05 | 00:20.9 | 19:00.0 | $0.0300$ | $0.0243$ | $0.0352$
$100.0$ | 2565:38.8 | 00:00.05 | 00:20.7 | 19:19.5 | $0.0303$ | $0.0249$ | $0.0352$

##  4 Conclusion

In this working paper, the potential of using PINNs to approximate forward
solutions of the Helmholtz equation is demonstrated for geometrically simple
2D and 3D problems (square or cubic computational domains) with homogeneous
Neumann and Dirichlet boundary conditions. Thereby, the readily available PINN
framework deepXDE is used to implement the neural networks and their
respective cost functions. Having very low training losses and excellent
agreement with the analytical solutions, the results exhibit a promising
ability of PINNs as forward solvers.

Regarding computational cost, Tab. 2 exhibits an interesting behavior: Using
the trained PINN as a forward-pass solver for the PDE, it is orders of
magnitude faster than using a FEM solver. This highlights the potential of
trained PINNs as surrogate models within optimization frameworks. However, it
has to be noted that a study regarding the generalization abilities of the
PINN is required prior to using the PINN in an optimization framework.

Future work may include a hyperparameter optimization of the PINNs in order to
achieve a smaller test loss. Additionally, it is subject to further research,
how well the PINN approach generalizes to more complex geometries to be
applicable in real-world problems, i.e. large complex non-cuboid shapes.
Nevertheless, the results presented in this working paper encourage future
work.

## Data and code availability

The data and code is available on reasonable request from the authors.

## References

  * [1] Johannes D Schmid, Philipp Bauerschmidt, Caglar Gurbuz, and Steffen Marburg.  Physics-informed neural networks for acoustic boundary admittance estimation.  2023, doi:10.2139/ssrn.4545338.  Preprint Available at SSRN 4545338.
  * [2] Florian Kraxberger, Eric Kurz, Werner Weselak, Gernot Kubin, Manfred Kaltenbacher, and Stefan Schoder.  A validated finite element model for room acoustic treatments with edge absorbers.  Acta Acustica, 7(48):1–19, October 2023, doi:10.1051/aacus/2023044.
  * [3] Cheol-Ho Jeong.  Room acoustic simulation and virtual reality-technological trends, challenges, and opportunities.  J. Swed. Acoust. Soc.(Ljudbladet), 1:27–30, 2022.
  * [4] Michael Vorländer.  Computer simulations in room acoustics: Concepts and uncertainties.  The Journal of the Acoustical Society of America, 133(3):1203–1213, 2013.
  * [5] Florian Kraxberger, Eric Kurz, Leon Merkel, Manfred Kaltenbacher, and Stefan Schoder.  Finite element simulation of edge absorbers for room acoustic applications.  In Fortschritte der Akustik — DAGA 2023, pages 1292–1295, Hamburg, March 2023. Deutsche Gesellschaft für Akustik.
  * [6] Florian Kraxberger, Eniz Museljic, Eric Kurz, Florian Toth, Manfred Kaltenbacher, and Stefan Schoder.  The nonlinear eigenfrequency problem of room acoustics with porous edge absorbers.  In Arianna Astolfi, Francesco Asdrubali, and Louena Shtrepi, editors, Proceedings of the 10th Convention of the European Acoustics Association: Forum Acusticum 2023, pages 6159–6166, Torino, September 2023. European Acoustics Association.
  * [7] Lu Lu, Xuhui Meng, Zhiping Mao, and George Em Karniadakis.  DeepXDE: A deep learning library for solving differential equations.  SIAM Review, 63(1):208–228, 2021, doi:10.1137/19M1274067.
  * [8] Stefan Schoder, Michael Weitz, Paul Maurerlehner, Alexander Hauser, Sebastian Falk, Stefan Kniesburges, Michael Döllinger, and Manfred Kaltenbacher.  Hybrid aeroacoustic approach for the efficient numerical simulation of human phonation.  The Journal of the Acoustical Society of America, 147(2):1179–1194, 2020.
  * [9] Stefan Schoder, Clemens Junger, and Manfred Kaltenbacher.  Computational aeroacoustics of the eaa benchmark case of an axial fan.  Acta Acustica, 4(5):22, 2020.
  * [10] Stefan Schoder, Klaus Roppert, Michael Weitz, Clemens Junger, and Manfred Kaltenbacher.  Aeroacoustic source term computation based on radial basis functions.  International journal for numerical methods in engineering, 121(9):2051–2067, 2020.
  * [11] Stefan Schoder, Manfred Kaltenbacher, Étienne Spieser, Hugo Vincent, Christophe Bogey, and Christophe Bailly.  Aeroacoustic wave equation based on pierce’s operator applied to the sound generated by a mixing layer.  In 28th AIAA/CEAS Aeroacoustics 2022 Conference, page 2896, 2022\.
  * [12] Lorenzo Tieghi, Stefan Becker, Alessandro Corsini, Giovanni Delibra, Stefan Schoder, and Felix Czwielong.  Machine-learning clustering methods applied to detection of noise sources in low-speed axial fan.  Journal of Engineering for Gas Turbines and Power, 145(3):031020, 2023.
  * [13] Stefan Schoder, Florian Kraxberger, Sebastian Falk, Andreas Wurzinger, Klaus Roppert, Stefan Kniesburges, Michael Döllinger, and Manfred Kaltenbacher.  Error detection and filtering of incompressible flow simulations for aeroacoustic predictions of human voice.  The Journal of the Acoustical Society of America, 152(3):1425–1436, 2022.
  * [14] Paul Maurerlehner, Stefan Schoder, Johannes Tieber, Clemens Freidhager, Helfried Steiner, Günter Brenn, Karl-Heinz Schäfer, Andreas Ennemoser, and Manfred Kaltenbacher.  Aeroacoustic formulations for confined flows based on incompressible flow data.  Acta Acustica, 6:45, 2022.
  * [15] Stefan Schoder, Étienne Spieser, Hugo Vincent, Christophe Bogey, and Christophe Bailly.  Acoustic modeling using the aeroacoustic wave equation based on pierce’s operator.  AIAA Journal, pages 1–10, 2023.
  * [16] Stefan Schoder and Andreas Wurzinger.  Dataset cylincf-01 creation pipeline: Circular cylinder in a cross flow, mach number 0.03 and reynolds number 200.  arXiv preprint arXiv:2303.05265, 2023.
  * [17] Andreas Wurzinger, Florian Kraxberger, Paul Maurerlehner, Bernhard Mayr-Mittermüller, Peter Rucz, Harald Sima, Manfred Kaltenbacher, and Stefan Schoder.  Experimental prediction method of free-field sound emissions using the boundary element method and laser scanning vibrometry.  Acoustics, 6(1):65–82, 2024.
  * [18] S. Schoder and K. Roppert.  openCFS: Open source finite element software for coupled field simulation – part acoustics, 2022.
  * [19] Maziar Raissi, Paris Perdikaris, and George E Karniadakis.  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, 2019.
  * [20] Stefan Schoder, Eniz Museljic, Florian Kraxberger, and Andreas Wurzinger.  Post-processing subsonic flows using physics-informed neural networks.  In Proceedings of 2023 AIAA Aviation Forum, pages 1–8, San Diego, June 2023. AIAA.
  * [21] Chao Song, Tariq Alkhalifah, and Umair Bin Waheed.  Solving the frequency-domain acoustic VTI wave equation using physics-informed neural networks.  Geophysical Journal International, 225(2):846–859, January 2021, doi:10.1093/gji/ggab010.
  * [22] Paul Escapil-Inchauspé and Gonzalo A. Ruz.  Hyper-parameter tuning of physics-informed neural networks: Application to helmholtz problems.  Neurocomputing, 561:126826, 2023, doi:10.1016/j.neucom.2023.126826.
  * [23] Rini J. Gladstone, Mohammad A. Nabian, and Hadi Meidani.  Fo-pinns: A first-order formulation for physics informed neural networks.  Arxiv Preprint, 2022, doi:10.48550/arxiv.2210.14320.
  * [24] Yanqi Wu, Hossein S. Aghamiry, Stephane Operto, and Jianwei Ma.  Helmholtz-equation solution in nonsmooth media by a physics-informed neural network incorporating quadratic terms and a perfectly matching layer condition.  Geophysics, 88(4):T185–T202, June 2023, doi:10.1190/geo2022-0479.1.
  * [25] Reza Akbarian Bafghi and Maziar Raissi.  Pinns-tf2: Fast and user-friendly physics-informed neural networks in tensorflow v2.  Arxiv Preprint, 2023, doi:10.48550/arXiv.2311.03626.
  * [26] Johannes Schmid, Philipp Bauerschmidt, Caglar Gurbuz, and Steffen Marburg.  Physics-informed neural networks for acoustic boundary admittance estimation.  Mechanical Systems and Signal Processing, 2023, doi:10.2139/ssrn.4545338.  Preprint.
  * [27] Frank Rosenblatt.  The perceptron: a probabilistic model for information storage and organization in the brain.  Psychological review, 65(6):386, 1958.
  * [28] Simon Haykin.  Neural networks: a comprehensive foundation.  Prentice Hall PTR, 1998.
  * [29] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [30] Dong C Liu and Jorge Nocedal.  On the limited memory bfgs method for large scale optimization.  Mathematical programming, 45(1-3):503–528, 1989.
  * [31] David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams.  Learning representations by back-propagating errors.  nature, 323(6088):533–536, 1986.
  * [32] Atilim Gunes Baydin, Barak A Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind.  Automatic differentiation in machine learning: a survey.  Journal of Marchine Learning Research, 18:1–43, 2018.


# Energy efficiency of slotted LoRaWAN communication with out-of-band
synchronization

Luca Beltramelli, Aamir Mahmood, Patrik Österberg, Mikael Gidlund, Paolo
Ferrari, and Emiliano Sisinni L. Beltramelli, A. Mahmood, P. Österberg, and M.
Gidlund are with the Department of Information Systems and Technology, Mid
Sweden University, 851 70 Sundsvall, Sweden, e-mail<EMAIL_ADDRESS>P. Ferrari and E. Sisinni are with the Department of Information Engineering,
University of Brescia, 25123 Brescia, Italy<EMAIL_ADDRESS>
###### Abstract

Although the idea of using wireless links for covering large areas is not new,
the advent of Low Power Wide Area Networks (LPWANs) has recently started
changing the game. Simple, robust, narrowband modulation schemes permit the
implementation of low-cost radio devices offering high receiver sensitivity,
thus improving the overall link budget. The several technologies belonging to
the LPWAN family, including the well-known LoRaWAN solution, provide a cost-
effective answer to many Internet-of-things (IoT) applications, requiring
wireless communication capable of supporting large networks of many devices
(e.g., smart metering). Generally, the adopted medium access control (MAC)
strategy is based on pure ALOHA, which, among other things, allows to minimize
the traffic overhead under constrained duty cycle limitations of the
unlicensed bands. Unfortunately, ALOHA suffers from poor scalability, rapidly
collapsing in dense networks. This work investigates the design of an improved
LoRaWAN MAC scheme based on slotted ALOHA. In particular, the required time
dissemination is provided by out-of-band communications leveraging on Radio
Data System (FM-RDS) broadcasting, which natively covers wide areas both
indoor and outdoor. An experimental setup based on low-cost hardware is used
to characterize the obtainable synchronization performance and derive a timing
error model. Consequently, improvements in success probability and energy
efficiency have been validated by means of simulations in very large networks
with up to 10000 nodes. It is shown that the advantage of the proposed scheme
over conventional LoRaWAN communication is up to 100% when short update time
and large payload are required. Similar results are obtained regarding the
energy efficiency improvement, that is close to 100% for relatively short
transmission intervals and long message duration; however, due to the
additional overhead for listening the time dissemination messages, efficiency
gain can be negative for very short duration of message fastly repeating.

###### Index Terms:

LoRaWAN, LPWAN, slotted ALOHA, time synchronization, FM-RDS.

##  I Introduction

The medium access control (MAC) mechanism is the most relevant part of the
data link layer; placed just above the physical layer, it controls how the
common network resources are shared. For this reason, it is one of the main
control knobs of any communication solution, which must fully exploit the
capabilities of the underlying physical layer without impairing application
requirements.

When typical (wireless) Internet-of-things (IoT) applications are considered
[1], (i.e., for transmitting data from many simple and autonomous devices
deployed in the field towards the Cloud), the emphasis is on simplicity and
security, while other aspects are somewhat overshadowed. Distributed MAC
approaches are generally preferred since they minimize the traffic overhead,
avoiding the need for forwarding control information towards a single node. In
particular, for non-critical IoT applications, random access paradigms are
often used. Although it is an inefficient mechanism, pure ALOHA is still
widely used due to its many advantages; packets can have variable size, nodes
can start transmission at any time, and time synchronization is not required.
For all these reasons, pure ALOHA has been selected as the MAC protocol for
most (if not all) Low Power Wide Area Networks (LPWANs), including the broadly
adopted LoRaWAN [2]. However, the main limitation of ALOHA is the large number
of collisions that can occur when the network consists of a large number of
devices. Better performance could be obtained using the ”listen-before-talk”
(LBT) approach, as in carrier sense multiple access with collision avoidance
(CSMA/CA). However, the effects of overhearing, hidden nodes, and long ongoing
transmission detection can negate the theoretical performance improvement for
low data rate communication over wide area.

The transmission vulnerability time of ALOHA can be halved, assuming all the
devices are relatively synchronized, by aligning the transmissions to slot
boundaries as in slotted ALOHA (S-ALOHA) [3]. Such an approach could be
particularly interesting for applications requiring periodic sampling of the
devices, which must be somehow time synchronized [4], as common in most IoT-
like monitoring applications such as smart metering [5]. Thus, the time-
dissemination strategy required by the application layer (which must schedule
the sampling activity), can also be used by the lower layers of the
communication protocol stack. As a matter of fact, some researchers have
already proposed the use of simple, in-band synchronization mechanisms for
enhancing LoRaWAN communication [6, 7, 8, 9].

However, the applications of interest mostly rely on uplink messages (from the
field toward a data sink) and should aim at limiting the traffic overhead,
especially in downlink (e.g., for exchanging time-related information), in
order to minimize the risk of collisions and increase the scalability. As a
consequence, the overhead imposed by in-band synchronization can result into
intolerable reduction of the communication opportunities for the devices. On
the contrary, out-of-band mechanisms could provide a more energy-efficient and
effective (i.e., offering high capacity/device density) approach under the low
duty-cycle restrictions imposed on LoRaWAN.

This work aims to propose and evaluate the performance of an S-ALOHA scheme
for LoRaWAN using an out-of-band synchronization technology. In particular,
the main contributions are:

  * •

analysis of time-dissemination technologies able to cover wide area including
indoor and outdoor scenarios with focus on FM-Radio Data System (FM-RDS),

  * •

theoretical modeling and experimental evaluation of the timing errors in a
real-world implementation,

  * •

comparative analysis of simulation results between LoRaWAN and the proposed
communication mechanism in terms of transmission success probability and
energy efficiency for large-scale deployment in an urban scenario.

The rest of the paper is organized as follows. The details of LoRaWAN and
time-dissemination technologies are presented in Sec. II. The proposed S-ALOHA
approach is described in Sec. III, while timing errors are analyzed in Sec.
IV. Experiments are detailed in Sec. V, while simulation results are reported
in Sec. VI. Finally, conclusions are drawn in Sec. VII.

##  II LoRaWAN and the out-of-band synchronization

As stated earlier, LPWANs are developed to provide multi-km network coverage
for wireless IoT applications. They typically employ single-hop connectivity
and are designed in order to minimize power consumption at the expense of data
rate (of few kbps) and of latency (in the order of seconds). From the
application point of view, the main target of LPWANs is to support delay-
tolerant massive machine-type communications (mMTC). Most mMTC applications
have some common traffic characteristics; in particular, each device generates
small and possibly infrequent messages, imposing a much higher traffic on
uplink [10]. On the contrary, with the exception for adaptive communication
parameters update, downlink is seldomly used. Examples of such uplink-centric
applications, implementing ”fire-and-forget” communication strategy, include
smart metering, fleet and asset tracking, vending machines management.
Currently, the term LPWAN includes many different technologies, standardized
or proprietary, operating in licensed or licence-free bands. The solutions
that can boast the widest adoption and provide similar features are arguably
NB-IoT, SigFox, and LoRa/LoRaWAN. Their most important characteristics are
summarized in Table I. For a comprehensive comparison of these technologies,
we refer the reader to [11]. For the sake of completeness, note that 3GPP also
supports LTE Cat-M, which is a solution fully compliant with LTE that
leverages on reduced channel bandwidth to offer raw throughput in the order of
few Mbps.

TABLE I: Salient Features of different LPWAN Solutions

Features |  Sigfox |  LoRaWAN |  NB-IoT (NB2)
---|---|---|---
Modulation |  BPSK |  CSS |  QPSK
Spectrum |  Unlicensed ISM bands |  Unlicensed ISM bands |  Licensed LTE bands
Bandwidth |  100 Hz |  250/125 kHz |  180 kHz
Max data rate |  100 bps |  50 kbps |  127 kbps (DL) / 159 kbps (UL)
Bidirectional |  Limited (Half-duplex) |  Yes (Half-duplex) |  Yes (Half-duplex)
Max messages/day |  140 (UL), 4 (DL) |  Duty-cycle constrained |  Unlimited
Max payload length |  12 bytes (UL), 8 bytes (DL) |  242 bytes |  2536 bit
Interference immunity |  Very high |  Very high |  Low
Authentication & encryption |  Yes (Enc. is optional, AES128) |  Yes (AES128) |  Yes (LTE enc.)
Adaptive data rate |  No |  Yes |  No
Handover |  End-devices do not join BS |  End-devices do not join to BS |  End-devices join single BS
Private nwk |  No |  Yes |  No
Standardization |  Sigfox / ETSI |  LoRa-Alliance |  3GPP

In the following, a brief overview of LoRaWAN is provided, followed by a
discussion on viable out-of-band time-dissemination technologies offering
wide-area coverage.

###  II-A LoRaWAN: A brief overview

LoRaWAN is a member of the LPWANs family, which includes all the communication
solutions aimed at transferring limited amount of data over a wide area with
relaxed time constraints. LPWANs generally operate in the sub-GHz unlicensed
bands (that provides better propagation compared to the crowded $2.4$ GHz
band) and exploit the reduced channel bandwidth, thus trading smaller data
rates for increased receiver sensitivity. LoRaWAN is based on a proprietary
radio technology (named LoRa), leveraging on an efficient chirp spread
spectrum modulation. The poor bandwidth availability is overcome by quasi-
orthogonal virtual channels, offering different data rates. These virtual
channels depend on the spreading factor (SF): a tunable parameter, which
modifies the symbol duration and represents the number of bits coded per
symbol. In particular, given the SF, the chirp length is defined (affecting
the data rate as well), and $2^{\textrm{SF}}$ different frequency trajectories
exist. Additional interleaving and whitening strategies are applied to
increase interference robustness, and forward error correction with coding
rate $\textrm{CR}\in[4/5,..,4/8]$ is adopted as well. LoRa offers so many
advantages that it is the foundation of many proprietary solutions (e.g., [12,
13]), other than LoRaWAN, which provides a standardized protocol stack.

The LoRaWAN channel access mechanism is pure ALOHA, chosen in order to
minimize stack complexity and traffic overhead for communication management.
However, it is well-known that ALOHA can be a significant bottleneck to the
LoRaWAN scalability in dense scenarios. For this reason, in the literature
some advanced random and scheduled access solutions have been discussed [9, 6,
14, 15]. The possible use of LBT, referred to as channel activity detection
(CAD), is rarely adopted due to the limits mentioned earlier. Scheduling of
transmissions reduces message collisions, however, it requires additional
overhead for providing time synchronization across the whole network.

From the topology point of view, LoRaWAN implements a star-of-stars
architecture, in which GWs forward every incoming wireless message to the
backhaul network and vice versa. The GWs do not manipulate user data, and the
node identification and authentication are managed by the Network Servers
(NS), whereas user data are managed by the Application Server (AS). LoRa
devices offer different features, depending on their class; Class-A permits
event-based uplink followed by an (optional) downlink, Class-B offers
synchronized downlink, lastly, Class-C implements continuous listening. For
the sake of clarity, in the rest of the paper, the term LoRaWAN refers to the
complete communication solution, whereas LoRa to the physical layer.

###  II-B Large areas out-of-band time-dissemination technologies

A node executes a synchronization algorithm to align the time of its local
clock to the clocks of other nodes or to a global clock. In particular,
absolute synchronization requires the distributed clocks to refer to an actual
global time standard (e.g., the Universal Coordinated Time UTC), whereas
relative synchronization requires time alignment within the local system,
involving only local nodes. In most practical cases, absolute synchronization
is not a mandatory requirement.

Synchronization implies transmission of time-related information and, when it
occurs over wireless networks, must respect many constraints dictated by
limited resource availability (energy, bandwidth, computation capability),
network topology, and unreliable links. The cost factor is also important when
the application scales. In the case of LPWANs, offering poor throughput and
reduced number of physical channels, the constraints become more stringent. As
a result, the overhead of two-way communication for in-band synchronization
mechanisms can unacceptably reduce communication opportunities (e.g., due to
limitations on the duty-cycle) [9] especially for the aforementioned ”fire-
and-forget”, uplink-centric, monitoring applications. For this reason, this
work evaluates the out-of-band approach (sometimes referred to as hardware-
assisted synchronization [16]), where a complementary communication system is
used (possibly wireless, as for wake-up receivers [17]). The criteria for the
selection of suitable synchronization technology for LPWANs are: a) accurate
event timestamping for minimizing the local clock offset with respect to the
chosen time reference (in other words, the complementary system must minimize
the uncertainty of the time of arrival estimation, b) wide-area coverage
complemented by indoor/outdoor capability, and c) minimal energy consumption
overhead.

In this work, we consider the following out-of-band technologies: global
navigation satellite system (GNSS), radio controlled clocks (RCCs), and the
Radio Data System (FM-RDS).

####  II-B1 GNSS - GPS

The Global Positioning System (GPS) is the most diffused GNSS, which relies on
strict time synchronization for evaluating pseudoranges in trilateration.
Nowadays, GPS receivers can provide a 1-PPS time reference signal while
offering excellent short-term stability and sub-microsecond synchronization
accuracy [18]. Unfortunately, the cost of GPS receivers may be too high for
several IoT applications, and their power consumption may not be compatible
with battery-operated devices [19]. Moreover, GPS receivers require line-of-
sight with satellites and thus cannot be used indoor.

####  II-B2 RCC - DCF77

The long-wave RCCs operate in the $40$ kHz to $80$ kHz frequency bands. The
use of low-frequency signaling offers wide-area coverage, lower power than
typical radio frequency signals, and better indoor reception. The DCF77
transmitter is operated by Physikalisch Technische Bundesanstalt in
Frankfurt/Main, which covers all the central Europe. In the DCF77 signal, the
beginning of each second can be easily recognized, providing a time accuracy
in the order of $0.1$ s. Using an additional phase modulated signal
superimposed on the amplitude modulation (AM) signal and using correlation
receiver, the synchronization accuracy can be improved to a few hundreds of
microseconds [20].

####  II-B3 FM-RDS

RDS is a supplementary digital data service that is superimposed on the
regular transmission of an FM broadcast station. RDS is intended to provide
additional, real-time services including, but not limited to, road traffic
information. The digital data are transmitted via double-sideband AM on a
suppressed $57$ kHz carrier, which is the third harmonic of the pilot tone for
FM stereo transmission. Data are biphase coded, and the data rate is $1187.5$
bps. Data are hierarchically arranged into 16-bit words, 26-bit blocks, and
104-bit groups. Depending on the data transmitted, FM-RDS defines multiple
types of groups, each identified by a 4-bit group type code. A clock time and
date (CT) group exists, purposely designed to transmit time information for
receiver synchronization. A CT-group is transmitted once per minute, but the
provided accuracy is poor, in the order of $100$ ms [21]. Unfortunately, the
CT may be absent in some RDS data stream, since it is not mandatory and some
broadcasting stations do not use it. Moreover, using FM-RDS for absolute
synchronization is unreliable, since broadcasting stations are not required to
be UTC synchronized. Subsequently, when group transmissions are considered,
relative time synchronization can be achieved with an error in the order of
hundreds of microseconds [22], however, offset with UTC may be unknown.
Finally, it has to be stressed that low-power consumption FM-RDS receivers can
be realized [23].

##  III The proposed communication mechanism

In this work, the use of slotted ALOHA is suggested as a viable approach in
dense LoRaWAN networks, under the name of S-LoRa protocol, to improve the
overall energy efficiency and scalability. Within the scope of this work, we
consider a monitoring application scenario with uplink only traffic. As stated
before, in such a case the need for downlink frames for transferring time-
related information can lead to excessive usage of communication
opportunities, already bounded by severe constraints in terms of available
bandwidth and duty-cycle limits. For this reason, differently from other works
in the literature, out-of-band synchronization is used. A common sense of time
among nodes is often required by the application to manage events of interests
(e.g., in monitoring applications, to periodically sample the quantity of
interest and schedule the transmission of the acquired values). In this work,
a common sense of time is added to the original LoRaWAN stack, and used to
define the timeslot boundaries within a contention window interval $T_{CW}$.
In particular, relative time synchronization (in which a synchronization event
simultaneously received by all the nodes is tracked) is considered. In detail,
the synchronization procedure requires:

  * •

two successive synchronization events, separated by a known interval, to
estimate the local clock rate drift with respect to the common reference;

  * •

a single synchronization event for correcting the clock phase offset error.

All the technologies discussed in Sec. II-B can be adopted; however, in this
work, we focus on the FM-RDS solution that, although not offering the best
accuracy, provides indoor coverage and ensures reliability due to the many FM
stations broadcasting in a single region.

As previously discussed in Sec. II-B3, CT-groups transmitted by an FM
broadcaster contain the time information, but the accuracy is poor.
Consequently, we propose to use the FM-RDS CT-group detection as a common
reference for the relative synchronization of the LoRa devices. The timestamp
contained in each CT-group is only used during the initialization procedure
(i.e., when a node joins the network). The only requirement is that all the
nodes must simultaneously recognize the same CT-group and accurately timestamp
it with the local clock. Fig.1 shows a flowchart of the proposed transmission
scheme, including the main steps required for a device to synchronize and
transmit a message to the gateway, compared to a regular event-based LoRaWAN
transmission.

Figure 1: Flowchart of: a) regular LoRaWAN transmissions, b) the proposed
S-LoRa approach based on out-of-band time dissemination.

To start the synchronization process, a device wakes up its FM-RDS receiver in
the proximity of the beginning of a contention window, compensating for the
accumulated clock drift, and waits for a CT-group. After a
$T_{\textrm{SYNC}}$, denoting the CT-group periodicity nominally equivalent to
one minute, the node receives a second CT-group, which allows estimating the
clock drift. The second CT-group also marks the beginning of the new
contention window. Once the synchronization is completed, the device waits for
a randomly chosen slot to transmit its message to the GW. The timeslot
duration (see Fig. 2) is designed to contain a message transmission, lasting a
time on air (ToA), plus a guard interval ($T_{g}$), i.e.,
$T_{\textrm{slot}}=ToA+T_{g}$, for avoiding the collisions between
transmissions in different slots. Let $T_{TX}=k\cdot T_{\textrm{SYNC}}$, for
some value $k\in\mathbb{Z}^{+}$, be the time interval available for
transmitting, the contention window has a duration $T_{CW}=T_{TX}-\Delta$,
where, $\Delta$ takes into account the non integer number of
$T_{\textrm{SYNC}}$ intervals (nominally lasting for one minute) encompassing
the actual contention window, and $T_{TX}$ is sized accordingly to the target
refresh period. Hence, the number of available slots is
$M=\frac{T_{CW}}{T_{\textrm{slot}}}=\frac{T_{TX}-\Delta}{ToA+T_{g}}$.

A message is successfully transmitted to the GW if there are no intra- and
inter-slot collisions. The former occurs when more than one device chooses the
same slot for transmission; the latter occurs when messages transmitted in
adjacent timeslots partially overlap due to synchronization errors and
transmission time uncertainties. Although not within the scope of this paper,
it is possible to derive the optimal guard time $T_{g}$, maximizing the
overall throughput.

Finally, it is worth noticing that, despite the MAC mechanism of the devices
has been modified, no LoRaWAN rule is violated and the message payload is
unaltered. Accordingly, the proposed approach does not require changes in the
legacy LoRaWAN backend infrastructure, ensuring low installation and
maintenance costs.

Figure 2: Structure of the timeslots for the proposed S-LoRa solution.
Synchronization process requires the reception of two CT-groups by the FM-RDS
radio.

##  IV Timing errors

For any synchronization mechanism, timing errors arise from the propagation
delay of the messages, and the non-ideal behavior of the transceivers and the
local clocks. In this section, the factors contributing to the uncertainty of
the transmission time of a device in the proposed communication mechanism are
discussed and characterized.

###  IV-A Transmission timing error

A transceiver is responsible for transmitting the LoRa frame once the timeslot
deadline has occurred. Unfortunately, traveling across the lower layers of the
node’s protocol stack adds delays, which cannot be easily entirely removed
even by proper calibration. Typically, such delays depend on the
software/hardware being used, and can be different from manufacturer to
manufacturer. Therefore, the Type A uncertainty $u_{TX}$ must be taken into
account, which we evaluated experimentally (see Sec. V-A).

###  IV-B Message propagation delay

Two messages transmitted precisely at the start of their respective timeslots
can still collide at the receiver due to the propagation delay. In the case of
line-of-sight conditions, the propagation delay increases linearly with the
distance of a device from the gateway. Assuming that the devices are uniformly
distributed in an annular shaped region, with radii $R_{L}$ and $R_{l}$, the
probability distribution $f_{PD}(x)$ of the propagation delay for a device is
given by

$\displaystyle f_{PD}\left(x\right)=$
$\displaystyle\begin{cases}\frac{2xv^{2}}{R_{L}^{2}-R_{l}^{2}}&\text{if
}\frac{R_{l}}{v}\leq x\leq\frac{R_{L}}{v}\text{,}\\\
0&\text{otherwise}\text{,}\end{cases}$ (1)

where $v$ is the speed of light. The distribution $f_{PD}(x)$ is a triangular
probability density function with mean $\mu_{PD}$ and standard deviation
$\sigma_{PD}$

$\displaystyle\mu_{PD}=\frac{2\left(R_{L}^{2}+R_{l}^{2}+R_{L}R_{l}\right)}{3\left(R_{L}+R_{l}\right)v}$
(2)
$\displaystyle\sigma_{PD}=\frac{\left(R_{L}-R_{l}\right)\sqrt{R_{L}^{2}+R_{l}^{2}+4R_{L}R_{l}}}{3\sqrt{2}v\left(R_{L}+R_{l}\right)}.$
(3)

For devices distributed on a disk-shaped region of radius $R$, we have
$R_{L}=R$ and $R_{l}=0$, giving $\mu_{PD}=\frac{2}{3}\frac{R}{v}$ and
$\sigma_{PD}=\frac{R}{3\sqrt{2}v}$. The Type B uncertainty $u_{PD}$ can be
estimated as

$u_{PD}=\sqrt{\mu_{PD}^{2}+\sigma_{PD}^{2}},$ (4)

the average value $\mu_{PD}$ is taken into account since it is not a-priori
known.

###  IV-C Local clock error

The natural or true clock is generally denoted by $t$ and runs at a rate of 1
second per second. The origin $t=0$ occurs at some arbitrary instant.
Unfortunately, real-world implementations of local clocks are imperfect, hence
their internal value is represented by the quantity $C(t)$ when evaluated at
the true instant $t$. Typically, a device clock is implemented through a timer
peripheral in a microprocessor-based system. Assuming the timer is running at
the nominal rate $F_{0}$, the timer period $T_{0}=1/F_{0}$ is the resolution
of $C(t)$, i.e., the smallest unit by which it is updated. Accordingly, when a
first-order clock model is adopted (sometimes referred to as the SKM model
[24]), $C(t)$ is described as

$C(t)\approx(\alpha+\beta t)T_{0},$ (5)

where the quantity $\alpha\cdot T_{0}$ is the time difference between local
clock time and true time at the origin, and $\beta\cdot T_{0}$ is the real
local clock rate. The clock drift $\gamma$ and offset $\theta(t)$ can be,
respectively, defined as

$\gamma=\frac{dC(t)}{dt}-F_{0}T_{0}=\frac{dC(t)}{dt}-1=(\beta-F_{0})T_{0},$
(6) $\theta(t)=C(t)-t.$ (7)

For a perfectly syntonized clock, the drift (or skew) is $\gamma=0$ and
$\beta=F_{0}$; for a perfectly synchronized clock the offset is $\theta(t)=0$
and $\alpha=0$.

As a consequence, from the node point of view, the start of the next
transmission (not necessarily the start of the next timeslot), $t_{d}$ in the
true time reference, is seen as the local clock value $C(t_{d})$. If a
synchronization procedure is carried out by the node, it is possible to assume
that the timer is restarted and the actual clock value is represented by

$V(t_{d})=C(t_{d})-C(t_{0})=\beta(t_{d}-t_{0})T_{0}\text{,}$ (8)

where $t_{0}$ is the true clock value at the origin of the current
synchronization interval.

Thus, the timing error of a node relative to the start of a message
transmission is the result of the three different contributing factors: the
uncertainty in the syntonization procedure ($u_{\beta}$) causing the
transmitter to drift from the nominal rate $F_{0}$; the uncertainty
$u_{t_{0}}$ in the $t_{0}$ instant (e.g., due to the detection mechanism of
the synchronization event); and the uncertainty $u_{t_{d}}$ in the
transmission time $t_{d}$. The overall uncertainty $u_{V}$ can be obtained
from (8) by propagating these contributions, resulting in

$u_{V}=\sqrt{\left(\frac{\partial{V}}{\partial{\beta}}\right)^{2}\\!\cdot
u_{\beta}^{2}+\left(\frac{\partial{V}}{\partial{t_{d}}}\right)^{2}\\!\cdot
u_{t_{d}}^{2}+\left(\frac{\partial{V}}{\partial{t_{0}}}\right)^{2}\\!\cdot
u_{t_{0}}^{2}}\text{,}$ (9)

which is finally evaluated as

$u_{V}=T_{0}\sqrt{\left(t_{d}-t_{0}\right)^{2}\cdot
u_{\beta}^{2}+\beta^{2}\cdot u_{t_{d}}^{2}+\beta^{2}\cdot
u_{t_{0}}^{2}}\text{.}$ (10)

It is possible to further elaborate the uncertainty contributions using Fig.
3, in which it is assumed that the devices can latch the synchronization
events (e.g., the reception of a new CT-group) and transmit a new message only
on the rising edge of the local clock $V(t)$ running at $F_{0}$. Accordingly,
the actual deadline $\hat{t_{d}}$ is affected by the time quantization error
uniformly distributed between 0 and $T_{0}$, resulting in the uncertainty
$u_{t_{d}}$ as

$u_{t_{d}}=\frac{T_{0}}{\sqrt{3}}.$ (11)

The actual reference time instant $\hat{t_{0}}$ is affected by two errors, i)
the time quantization, as in (11), and ii) the timing error in the detection
of the synchronization event. These errors are characterized by the
uncertainties $u_{t_{0}q}$ and $u_{t_{0}s}$, respectively, leading to

$u_{t_{0}}=\sqrt{u_{t_{0}q}^{2}+u_{t_{0}s}^{2}}=\sqrt{\frac{T_{0}}{3}^{2}+u_{t_{0}s}^{2}}.$
(12)

Figure 3: Timing errors diagram highlighting the contributions of $u_{t_{0}s}$
and $u_{t_{0}q}$ to the clock synchronization uncertainty, and $u_{t_{0}d}$ to
the message transmission uncertainty.

Finally, the syntonization mechanism must be analyzed for estimating the
uncertainty $u_{\beta}$.

The $T_{\textrm{SYNC}}$ comprises the interval between two CT-group
transmission, i.e., $T_{\textrm{SYNC}}=60$ s. Consequently, the rate $\beta$
can be estimated as

$\hat{\beta}=\frac{\Delta_{\textrm{COUNTER}}}{T_{\textrm{SYNC}}},$ (13)

where $\Delta_{\textrm{COUNTER}}$ is a value representing the number of clock
periods $T_{0}$ counted in $T_{\textrm{SYNC}}$. Using (13) to propagate the
uncertainty of $T_{\textrm{SYNC}}$ (i.e., $u_{T_{\textrm{SYNC}}}$) leads to

$u_{\beta}=\frac{\Delta_{\textrm{COUNTER}}}{T_{\textrm{SYNC}^{2}}}\cdot
u_{T_{\textrm{SYNC}}}=\frac{\beta}{T_{\textrm{SYNC}}}\cdot
u_{T_{\textrm{SYNC}}}.$ (14)

The estimation $\hat{T}_{\textrm{SYNC}}$ of the synchronization interval is
affected by the time quantization, as well as the error in the detection of
the reference time signal. Differently from $\hat{t_{0}}$, their effects occur
both at the start and at the end of the interval itself. Consequently, the
uncertainty is

${u_{T_{\textrm{SYNC}}}=\sqrt{4u_{t_{0}q}^{2}+4u_{t_{0}s}^{2}}=\sqrt{4\frac{T_{0}}{3}^{2}+4u_{t_{0}s}^{2}}.}$
(15)

When the worst case scenario is considered, $t_{d}-t_{0}\simeq(k\cdot
T_{\textrm{SYNC}})=T_{TX}$ for some $k$, so that we can simplify $u_{V}$ as

$\displaystyle u_{V}$ $\displaystyle=\beta T_{0}\sqrt{k^{2}\cdot
u_{T_{\textrm{SYNC}}}^{2}+u_{t_{d}}^{2}+u_{t_{0}}^{2}}$ (16)
$\displaystyle=(1+\gamma)\sqrt{(4k^{2}+2)\frac{T_{0}}{3}^{2}+(4k^{2}+1)u_{t_{0}s}^{2}}.$

Assuming $k>>1$, $T_{0}<<u_{t_{0}s}$ and $\gamma<<1$, as typically occurs when
synchronization via FM-RDS signal is considered, (17) reduces to

$\displaystyle{u_{V}=2k\cdot u_{t_{0}s}.}$ (17)

Although $T_{0}$ is generally known, $u_{t_{0}s}$ must be evaluated, e.g., by
performing experiments in agreement with Type A uncertainty estimation. Such
an approach is adopted in this work, as detailed in Sec. V.

###  IV-D Concluding remarks

All the previously addressed errors contribute to the overall timing error
introduced in this section. Without loss of generality, the errors can be
grouped into the overall timing uncertainty $u$ given by

$u=\sqrt{u_{TX}^{2}+u_{PD}^{2}+u_{V}^{2}}.$ (18)

The value of $u$ can be transformed into an equivalent standard deviation of a
probability density function describing the possible timing error population
and used in the simulations of the considered application scenario. It has to
be emphasized that the proposed approach relies on the simultaneous estimation
of the local clock rate by all the devices using the same synchronization
events, so that variations in $T_{\textrm{SYNC}}$ affect all devices, without
compromising the timeslot boundaries.

##  V Experimental results

The accurate characterization of IoT components is known to be a main concern
of current research activities, requiring specific instrumentation setup and
measurement methods [25]. In this section, we describe the experimental
validation of hardware-assisted synchronization required for implementing the
S-LoRa protocol. We purposely designed an experimental setup with twofold
objectives:

  * •

to evaluate the uncertainty in transmitting a LoRaWAN message in a well-
defined time instant as needed by a timeslot deadline,

  * •

to evaluate the uncertainty in retrieving time information from the FM-RDS
broadcast signal.

In the considered scenario, focused on local relative time synchronization,
the propagation delay is generally negligible, given that devices are a few
kilometers apart from the gateway and the $u_{PD}$ is in the order of a few
microseconds, no matter the adopted SF. Each network node uses daughter radio
cards, which are connected to a host system (LoRa via an SPI link, and FM via
I2C link), executing the upper layers of the protocol stack. The LoRa radio is
based on embedded SX1272111Datasheet available at www.semtech.com. See also
[26]. daughter board, which is capable of signaling the completed transmission
(TX_DONE) or reception (RX_DONE) of LoRa frames by a digital line on the board
connector. This line is used to characterize the LoRa transceiver without the
additional latency caused by the execution of the protocol stack on the host
system.

Regarding the FM-RDS decoder, the Si4703 from Silicon Labs222Datasheet
available at www.silabs.com is used, hosted in an FM Radio Tuner Evaluation
Breakout Board from Hitletgo. This device implements a complete solution for
decoding both the audio and RDS signals, thus including many unused blocks
(e.g., the driver for the headphones); accordingly, it is not possible to
evaluate precisely the power consumption. When the ”verbose” mode is active, a
hardware pulse is generated each time a new RDS group is received. This
hardware line represents the actual time reference with respect to the
broadcasting station. The host board is an STM32 Nucleo board, running the
upper layers of the protocol stack (other than analyzing the content of the
incoming RDS CT-group).

###  V-A Characterization of LoRa frame transmission latency

We used two nodes (A and B) to characterize the capability of the LoRa radio
to timely transmit a new frame. In Node A, the Nucleo board is configured to
start the transmission after the recognition of a rising edge on its input
line associated with the highest interrupt priority. The trigger signal was
provided by a function generator (Keysight 33220A). To ensure a deterministic
behavior and minimize the overall latency, the frame content was already
available in memory before the trigger occurred. The Node B is permanently
configured in the reception mode to detect the LoRa frame. A counter (Keysight
53230A, receiving the $10$ MHz reference signal from the function generator)
is used (see Fig. 4) for measuring the time $\Delta_{TX}$ elapsing between the
trigger signal edge and the TX_DONE edge in Node A (configuration-a), and
between the trigger signal edge and the RX_DONE edge in the Node B
(configuration-b). The $\Delta_{TX}$ mean value reflects message duration and
hardware fixed delay (that can be removed by proper calibration), while the
standard deviation is $\sigma_{\Delta_{TX},a}=1.4~{}\mu$s in configuration-a
and $\sigma_{\Delta_{TX},b}=2.8~{}\mu$s in configuration-b, when uncertainties
of both Node A and Node B are combined; such values stay the same for
different SFs. In this way, it has been possible to estimate the uncertainty
in transmitting a LoRa frame $u_{TX}=1.4~{}\mu$s. For the sake of
completeness, note that the resolution of instantaneous frequency jumps inside
a single LoRa symbol is $T_{jump}=T_{C}/2^{\textrm{SF}}$, which is about
$8~{}\mu$s for $\textrm{SF}=7$.

Figure 4: Testbed for evaluating LoRa frame transmission latency.

###  V-B Characterization of FM-RDS CT-group reception

Two (co-located) nodes are used to characterize the capability of the FM
receivers to detect an incoming CT-group. Each node follows the FM-RDS data
stream continuously looking for CT-groups. Once a new CT-group is detected, a
pulse is generated on an STM32 digital output line ($CT_{OUT1,2}$ in Fig. 5).
The digital lines are acquired by a GPS-synchronized time server T103 from
HEOL, which is able to timestamp the two line edges with sub-microsecond
accuracy with respect to the UTC time. Initially, the population
$T_{CT_{OUT1}}$ of $CT_{OUT1}$ timestamps is considered for evaluating the
jitter in the CT-group periodicity, i.e., for evaluating the uncertainty
$u_{CT}$ in the time interval (nominally equal to $T_{\textrm{SYNC}}$)
elapsing between two successive CT-groups. In particular, we consider the
standard deviation $\sigma_{CT_{OUT1}}=u_{CT}$. Recall that $u_{CT}$ does not
affect the proposed synchronization strategy. Subsequently, the time
difference of arrival of the pulse at the nodes’ digital output lines
$\Delta_{CT}=T_{CT_{OUT1}}-T_{CT_{OUT2}}$ is considered for taking into
account the uncertainties of both nodes, which ideally would be equal. Because
of systematic errors, summarized by the average value $\mu_{\Delta_{CT}}$, the
CT-group reception uncertainty is estimated as

Figure 5: Testbed for evaluating CT-group message latency.

$\displaystyle
u_{t_{0}s}=\sqrt{\frac{\mu_{\Delta_{CT}}^{2}+\sigma_{\Delta_{CT}}^{2}}{2}},$
(19)

where $\sigma_{\Delta_{CT}}$ is the standard deviation of the $\Delta_{CT}$
population.

The value of $u_{CT}$ depends primarily on the mechanism employed by the FM
broadcaster for the transmission of CT-groups (in agreement with FM-RDS
specifications) and for the synchronization with the UTC absolute time. A two-
day long measurement campaign was carried out, syntonizing both the Si4703
devices on three different Italian FM broadcasting radio-stations. The “Radio
Monte Carlo” station achieved a $u_{CT}=171$ ms, the “RTL102.5” station showed
a $u_{CT}=62$ ms, and the “Radio Freccia” obtained a $u_{CT}=114$ ms. Only,
the best result is considered for the following part of the paper, thus
$u_{CT}=62$ ms is used in the simulation of Sec. VI. Additionally, for
“RTL102.5”, $\mu_{\Delta_{CT}}=0.41$ ms and $\sigma_{\Delta_{CT}}=0.24$ ms, so
that according to our experiments $u_{t_{0}s}=0.34$ ms.

##  VI Simulator Results

To evaluate the performance of the proposed slotted communication based on
out-of-band synchronization for a wide area metering application, we resorted
to using simulations. Indeed, simulations permit to easily and effectively
evaluate networks with thousands of devices. The performance is characterized
in terms of transmission success probability and energy efficiency, the latter
being the primary concern for any IoT wireless application.

###  VI-A The simulator

In this work, to simulate the proposed communication mechanism, we modified
LoRaEnergySim [27] Python simulator to incorporate the experimental and
analytical model results presented in Sec. IV and Sec. V333Simulator code
available at https://github.com/Beltra90/LoRaMACSim. Although numerous
LoRa/LoRaWAN simulators have been designed and reported in the literature (for
detailed comparison see e.g., [28]), careful consideration of the power
capture effect in each simulator is required for evaluating any new proposals.
In general, these simulators are based on two alternative power capture
models, either cumulative or dominant interference. For scenarios with many
interfering devices, the dominant interference model adopted by LoRaEnergySim
provides an upper bound to the performance of LoRa compared to the cumulative
interference model proposed in [29]. Consequently, when it is used to evaluate
the performance improvements of the proposed slotted communication, which
intrinsically reduces the number of interfering signals compared to the non-
slotted LoRaWAN, a dominant interference model provides more conservative
results. The salient features of the interference model of LoRaEnergySim are
as follows: a) signals transmitted using different SFs are orthogonal, b)
messages can only be correctly received if they satisfy the minimum signal-to-
interference ratio (SIR) and signal-to-noise ratio (SNR) thresholds, c) when
computing the SIR, only the dominant interferer is considered, and d) messages
can potentially survive a partial overlap limited to the first few symbols of
the preamble.

To model the time uncertainty affecting a message transmission, we opted to
use the worst-case combined uncertainty $u$ determined in Sec. IV-C. Each
transmission is affected by an independent and identically distributed
stochastic offset. The standard deviation of the timing offset distribution
affecting a transmission is chosen to be equal to the overall timing
uncertainty found in Sec. IV-D, (18). This model of the timing errors via the
combined uncertainty allowed us to simulate scenarios with a large number of
devices, avoiding the computationally expensive task of individually modeling
the local clock of each device.

###  VI-B Simulated scenarios

We simulated a typical urban deployment scenario in which a single GW is
located at the center of a group of devices. The devices are randomly and
uniformly distributed over a circular-shaped region with a diameter of $3$ km.
The signal propagation in the urban scenario was modeled similarly to [29],
however, we opted to use an uncorrelated lognormal shadowing fading with
standard deviation of $7.8$ dB [30]. The parameters used in the simulation are
reported in Table II. A convenient, linearly increasing, guard time $k\cdot
T_{g0}$ (where $k$ is the number of $T_{\textrm{SYNC}}$ in an $T_{TX}$) was
selected to maximize the throughput given the synchronization uncertainty $u$.
We tested different configuration of device densities, SFs, payload sizes and
transmission intervals $T_{TX}$, always respecting the limitations imposed by
the LoRaWAN specifications. The traffic in the scenario consists of
unconfirmed uplink messages, i.e., from end devices to the GW. Each device
generates a message once every $T_{TX}$, and the cumulative number of messages
transmitted by all devices during each simulation run is $200000$.
Consequently, in each run, the network is simulated for a time duration
ranging from one to 90 hours. For each combination of SF, $T_{TX}$, and
payload length $L_{P}$, the results are averaged over ten simulation runs. The
complete list of parameters used in the simulation is available in Table II.

TABLE II: LoRa Parameters used for the performance evaluation

Parameter | Symbol | Value
---|---|---
Bandwidth and Channel | BW and CH |  $125$ kHz and $868.1$ MHz
Preamble |  |  $8$ bytes
Coding Rate | CR | $4/8$
Spreading Factor | SF | $\in\left\\{7,9,12\right\\}$
Payload @ SF7 | $L_{P7}$ |  $\in\left\\{10,51,221\right\\}$ bytes
Payload @ SF9 | $L_{P9}$ |  $\in\left\\{10,51,115\right\\}$ bytes
Payload @ SF12 | $L_{P12}$ |  $\in\left\\{10,51\right\\}$ bytes
Guard time | $T_{g0}$ |  $3$ ms
Noise Figure and PSD |  |  $6$ dB and $-174$ dBm/Hz
Shadow fading |  |  $7.8$ dB
GW Antenna Height |  |  $15$ m
Device floor no. |  | $\mathcal{U}(1,4)$
Internal walls no. |  | $\mathcal{U}(0,3)$
SNR Thresholds |  |  $\left[-6,-12,-20\right]$ dB
SIR Threshold |  |  $1$ dB
Transmitting Power |  |  $14$ dBm
CT-group duration | $T_{CT}$ |  $86.7$ ms
Radio switch on time | $T_{ON}$ |  $1$ ms

###  VI-C Results: success probability and energy efficiency

The propoposed S-LoRa has been compared against the pure ALOHA-based
transmission of a LoRaWAN Class-A device transmitting unconfirmed messages. We
selected the (transmission) success probability as first metric of interest,
calculated as the fraction of generated messages that are successfully
received by the GW. For the set of parameters under investigation, the success
probability of S-LoRa and LoRaWAN was always observed to be smaller than 0.5.
It is possible to achieve significantly higher values by decreasing the
average distance of the devices from the gateway. However, this type of
consideration is beyond the scope of this paper, which instead primarily
focuses on the relative comparison of the performance of LoRaWAN and S-LoRa.
For this reason, in Fig. 6 the relative success probability gain of S-LoRa
with respect to the reference LoRaWAN ALOHA-like implementation is shown, for
a network consisting of $N=5000$ (Fig. 6.a), and $N=10000$ (Fig. 6.b) devices.
Depending on the transmission interval, ToA, and SF considered in the
scenario, improvements close to 100% in the success ratio can be reached,
confirming the feasibility and advantages of the proposed approach. As the
scenario becomes increasingly demanding, i.e., short $T_{TX}$ and large
$L_{P}$, the gain in success probability of S-LoRa over LoRaWAN becomes more
significant.

(a) $N=5000$ nodes in the network  (b) $N=10000$ nodes in the network

Figure 6: Percentage relative success probability gain of S-LoRa over LoRaWAN.

The other performance metric considered by this study is the energy
efficiency. It is important to investigate the energy efficiency of S-LoRa to
determine if the increased energy consumption resulting from the
synchronization mechanism is adequately compensated by an improvement in the
communication performance metric (i.e., the transmission success probability).
In calculating the energy consumption and energy efficiency the following
assumptions have been made:

  * •

both regular LoRaWAN and S-LoRa require an always-on clock for scheduling
monitoring activity and subsequent message transmission; therefore, the power
consumption of the clock oscillator circuit was not taken into account in the
comparison,

  * •

each time a device wants to synchronize, two successive CT-groups must be
received. The FM-RDS receiver remains active for $T_{RX1,\textrm{FM-
RDS}}=T_{ON}+T_{CT}+T_{TX}\cdot\gamma+\frac{3}{\sqrt{2}}\cdot u_{CT}$ in order
to receive the first synchronization event, where $T_{ON}$ and $T_{CT}$ are
defined in Table II. It is assumed that when the FM-RDS receiver is activated
after a long time the device will account for the uncertainty of its clock
drift. To receive the second CT-group, the FM-RDS remains active for
$T_{RX2,\textrm{FM-RDS}}=T_{ON}+T_{CT}+3\cdot u_{CT}$, since in this case the
device has adjusted its clock after the reception of the previous CT-group.

(a) $N=5000$ nodes in the network  (b) $N=10000$ nodes in the network

Figure 7: Percentage relative energy efficiency gain of S-LoRa over LoRaWAN.

The LoRa Model Calculator Tool from Semtech [26] is used to derive the
transmission and sleep currents for a SX1272/73 LoRa transceiver. Since
downlink traffic is not considered (only unconfirmed uplink messages are
sent), it is assumed that the LoRa transceiver is either in transmitting or
sleeping state and never in receiving state. The current in transmission state
is $I_{TX,\textrm{LoRa}}=44$ mA, whereas the current during sleep is
$I_{\textrm{sleep},\textrm{LoRa}}=100$ nA. Accordingly, assuming a supply
voltage $V_{\text{dd}}=3.3$ V, the power consumption during message
transmission is $P_{TX,\textrm{LoRa}}=145.2$ mW, which decreases to
$P_{\textrm{sleep},\textrm{LoRa}}=0.33~{}\mu$W in sleep state. Unfortunately,
FM-RDS receivers available on the market generally embed many functionalities
other than RDS stream demodulation, thus the current consumption reported in
many data sheets is significantly larger than the current strictly required
for the demodulation of the RDS stream. We used [23] as reference for the FM-
RDS receiver’s current consumption while active, where it is reported that the
absorbed current is $I_{RX,\textrm{FM-RDS}}=1.2$ mA with a supply voltage
$V_{\text{dd}}=3.3$ V, so that the power consumption during message
transmission is $P_{RX,\textrm{FM-RDS}}=3.96$ mW. When the FM-RDS receiver is
not listening for CT-groups, it is assumed to be in idle mode; due to
unavailability of current consumption in idle state in the literature, we
consider instead a state-of-the-art device as the CC1125 Narrowband
Transceiver [31], with a current consumption of $I_{\textrm{idle},\textrm{FM-
RDS}}=120$ nA and a supply voltage of $V_{\text{dd}}=3.3$ V, so that the power
consumption in the idle state is $P_{\textrm{sleep},\textrm{FM-
RDS}}=0.36~{}\mu$W. We obtained the switching times between different radio
states of the FM-RDS receiver from the same literature source, which we
aggregated and upper bounded by the quantity $T_{ON}$.

Relying on the previously described real-world hardware-related figures of
merit, the energy efficiency of a device is calculated as the ratio of the
overall number of bits successfully transmitted over the total energy
consumption of the device during the entire simulation. In Fig. 7, the
relative energy efficiency gain with respect to the reference LoRaWAN
implementation is shown, when the considered network consists of $N=5000$
(Fig. 7.a), and $N=10000$ (Fig. 7.b) devices. Like for the success probability
gain, we can observe that the gain in energy efficiency of S-LoRa over LoRaWAN
becomes more significant for short transmission intervals and large payload
sizes. We can also observe that although S-LoRa offers a strictly positive
gain in the success probability over LoRaWAN, the energy efficiency gain can
be negative (i.e., using S-LoRa can lead to a loss in energy efficiency).
Despite the proposed approach generally outperforms the ALOHA-like approach of
LoRaWAN in term of transmission success probability, from Fig. 7 we can
conclude that for applications requiring sporadic transmission of small
payloads, S-LoRa can lead to a reduction in energy efficiency, even in the
case of very crowded LoRa networks. This is especially the case for the
devices that, thanks to an advantageous link budget, can transmit at higher
data rates using the smaller SFs.

##  VII Conclusion

Currently, LoRaWAN represents one of the most widely adopted LPWAN technology,
capable of providing the long-range wireless connectivity required by many IoT
applications. The MAC layer in LoRaWAN, like that of many other LPWANs, is
based on pure ALOHA for minimizing the traffic overhead and for coping with
duty-cycle limitation. To improve the scalabilty of a single LoRaWAN gateway,
this work proposed the use of an improved LoRaWAN MAC scheme based on slotted
ALOHA (S-LoRa). Compared to related work in the literature, the main
contribution is the use of an out-of-band synchronization based on FM-RDS
broadcasting, a solution that natively covers wide areas both indoor and
outdoor. The paper provides an analytical model and an experimental
characterization of the synchronization performance of the proposed approach.
An extensive set of simulations, parameterized accordingly to a smart city
scenario (e.g., smart metering), was carried out to study the proposed
communication mechanism. In a dense network with 10000 devices, the proposed
scheme shows an improvement of up to 100% in terms of success probability
compared to the standard LoRaWAN approach. In addition, in terms of energy
efficiency, although the proposed S-LoRa communication can have a negative
gain in some conditions, it remains effective for relatively short
transmission intervals and large message sizes. As a concluding remark, it has
to be highlighted that the proposed time slot arrangement (actually
constituting what is normally defined as a superframe), can be modified in
order to include both uplink and downlink periods. Future works are planned
for evaluating hybrid in-band and out-of-band sync strategies, depending on
the overhead in terms of lost communication opportunities and increased power
consumption for possible overhearing.

## References

  * [1] A. Kumar, M. Zhao, K. Wong, Y. L. Guan, and P. H. J. Chong, “A comprehensive study of IoT and WSN MAC protocols: Research issues, challenges and opportunities,” _IEEE Access_ , vol. 6, pp. 76 228–76 262, 2018.
  * [2] Q. Song, X. Lagrange, and L. Nuaymi, “Evaluation of macro diversity gain in Long Range ALOHA networks,” _IEEE Commun. Lett._ , vol. 21, no. 11, pp. 2472–2475, 2017.
  * [3] ——, “An analytical model for S-ALOHA performance evaluation in M2M networks,” in _IEEE ICC_ , 2017, pp. 1–7.
  * [4] M. Rizzi, P. Ferrari, A. Flammini, and E. Sisinni, “Evaluation of the IoT LoRaWAN solution for distributed measurement applications,” _IEEE Trans. Instrum. Meas._ , vol. 66, no. 12, pp. 3340–3349, 2017.
  * [5] J. Lloret, J. Tomas, A. Canovas, and L. Parra, “An integrated IoT architecture for smart metering,” _IEEE Commun. Mag._ , vol. 54, no. 12, pp. 50–57, 2016.
  * [6] B. Reynders, Q. Wang, P. Tuset-Peiro, X. Vilajosana, and S. Pollin, “Improving reliability and scalability of LoRaWANs through lightweight scheduling,” _IEEE Internet Things J._ , vol. 5, no. 3, pp. 1830–1842, 2018.
  * [7] T. Polonelli, D. Brunelli, A. Marzocchi, and L. Benini, “Slotted ALOHA on LoRaWAN-design, analysis, and deployment,” _Sensors_ , vol. 19, no. 4, p. 838, 2019.
  * [8] M. Rizzi, A. Depari, P. Ferrari, A. Flammini, S. Rinaldi, and E. Sisinni, “Synchronization uncertainty versus power efficiency in LoRaWAN networks,” _IEEE Trans. Instrum. Meas._ , vol. 68, no. 4, pp. 1101–1111, 2019.
  * [9] J. Haxhibeqiri, I. Moerman, and J. Hoebeke, “Low overhead scheduling of LoRa transmissions for improved scalability,” _IEEE Internet Things J._ , vol. 6, no. 2, pp. 3097–3109, 2018.
  * [10] H. Shariatmadari, R. Ratasuk, S. Iraji, A. Laya, T. Taleb, R. Jäntti, and A. Ghosh, “Machine-type communications: current status and future perspectives toward 5G systems,” _IEEE Commun. Mag._ , vol. 53, no. 9, pp. 10–17, 2015.
  * [11] K. Mekki, E. Bajic, F. Chaxel, and F. Meyer, “A comparative study of LPWAN technologies for large-scale IoT deployment,” _ICT Express_ , vol. 5, no. 1, pp. 1–7, 2019.
  * [12] S. Gao, F. Shang, and C. Du, “Design of multichannel and multihop low-power wide-area network for aircraft vibration monitoring,” _IEEE Trans. Instrum. Meas._ , vol. 68, no. 12, pp. 4887–4895, 2019.
  * [13] H. Lee and K. Ke, “Monitoring of large-area IoT sensors using a LoRa wireless mesh network system: Design and evaluation,” _IEEE Trans. Instrum. Meas._ , vol. 67, no. 9, pp. 2177–2187, 2018.
  * [14] L. Leonardi, F. Battaglia, G. Patti, and L. L. Bello, “Industrial LoRa: A novel medium access strategy for LoRa in industry 4.0 applications,” in _IEEE IECON_ , 2018, pp. 4141–4146.
  * [15] K. Q. Abdelfadeel, D. Zorbas, V. Cionca, and D. Pesch, “${FREE}$–fine-grained scheduling for reliable and energy-efficient data collection in LoRaWAN,” _IEEE Internet Things J._ , vol. 7, no. 1, pp. 669–683, 2020.
  * [16] W. Ikram, I. Stoianov, and N. F. Thornhill, “Towards a radio-controlled time synchronized wireless sensor network: A work in-progress paper,” in _IEEE ETFA_ , Sep. 2010, pp. 1–4.
  * [17] F. A. Aoudia, M. Gautier, M. Magno, M. L. Gentil, O. Berder, and L. Benini, “Long-short range communication network leveraging LoRa™ and wake-up receiver,” _Microprocessors and Microsystems_ , vol. 56, pp. 184 – 192, 2018\.
  * [18] H. Guo and P. Crossley, “Design of a time synchronization system based on GPS and IEEE 1588 for transmission substations,” _IEEE Trans. Power Del._ , vol. 32, no. 4, pp. 2091–2100, Aug 2017.
  * [19] L. A. Tawalbeh, A. Basalamah, R. Mehmood, and H. Tawalbeh, “Greener and smarter phones for future cities: Characterizing the impact of GPS signal strength on power consumption,” _IEEE Access_ , vol. 4, pp. 858–868, 2016.
  * [20] A. Depari, A. Flammini, M. Lavarini, and E. Sisinni, “Inexpensive SDR-based longwave radio controlled clock for time dissemination in industrial wireless sensor networks,” in _IEEE Int. Conf. on Automation Science and Engineering (CASE)_ , Aug 2015, pp. 125–130.
  * [21] D. Kopitz and B. Marks, _RDS: the Radio Data System_. Boston, MA: Artech House, 1999.
  * [22] L. Li, L. Sun, G. Xing, W. Huangfu, R. Zhou, and H. Zhu, “ROCS: Exploiting FM radio data system for clock calibration in sensor networks,” _IEEE Trans. Mobile Comput._ , vol. 14, no. 10, pp. 2130–2144, 2015.
  * [23] D. Ruffieux, T. Melly, V. Peiris, J. . Perotto, N. Raemy, and E. Le Roux, “A 1.2 mW RDS receiver for portable applications,” _IEEE J. Solid-State Circuits_ , vol. 39, no. 7, pp. 995–1005, July 2004.
  * [24] D. Veitch, S. Babu, and A. Pàsztor, “Robust synchronization of software clocks across the internet,” in _ACM IMC_ , New York, NY, USA, 2004, pp. 219–232.
  * [25] M. Zhong, Y. Yang, H. Yao, X. Fu, O. A. Dobre, and O. Postolache, “5G and IoT: Towards a new era of communications and measurements,” _IEEE Instrum. Meas. Mag._ , vol. 22, no. 6, pp. 18–26, 2019.
  * [26] _SX127X2/3/6/7/8: LoRa Modem Designer’s Guide_ , Semtech, 2020 (Accessed on August). [Online]. Available: https://lora-developers.semtech.com/library/product-documents/
  * [27] G. Callebaut, G. Ottoy, and L. van der Perre, “Cross-layer framework and optimization for efficient use of the energy budget of IoT nodes,” in _IEEE WCNC_ , 2019, pp. 1–6.
  * [28] J. Haxhibeqiri, E. De Poorter, I. Moerman, and J. Hoebeke, “A survey of LoRaWAN for IoT: From technology to application,” _Sensors_ , vol. 18, no. 11, p. 3995, 2018.
  * [29] D. Magrin, M. Centenaro, and L. Vangelista, “Performance evaluation of LoRa networks in a smart city scenario,” in _IEEE ICC_ , 2017, pp. 1–7.
  * [30] J. Petajajarvi, K. Mikhaylov, A. Roivainen, T. Hanninen, and M. Pettissalo, “On the coverage of LPWANs: range evaluation and channel attenuation model for LoRa technology,” in _14th International Conference on ITS Telecommunications (ITST)_ , 2015, pp. 55–59.
  * [31] _Ultra-high performance Sub-1 GHz wireless transceiver for narrowband systems_ , Texas Instruments, 2020 (Accessed on August). [Online]. Available: https://www.ti.com/document-viewer/CC1125/datasheet