# Physics-Informed Neural Network based inverse framework for time-fractional
differential equations for rheology

Sukirt Thakur  School of Mechanical Engineering,Purdue University, West
Lafayette, 47907, Indiana, USA Harsa Mitra School of Mechanical
Engineering,Purdue University, West Lafayette, 47907, Indiana, USA Arezoo M.
Ardekani School of Mechanical Engineering,Purdue University, West Lafayette,
47907, Indiana, USA

###### Abstract

Time-fractional differential equations offer a robust framework for capturing
intricate phenomena characterized by memory effects, particularly in fields
like biotransport and rheology. However, solving inverse problems involving
fractional derivatives presents notable challenges, including issues related
to stability and uniqueness. While Physics-Informed Neural Networks (PINNs)
have emerged as effective tools for solving inverse problems, most existing
PINN frameworks primarily focus on integer-ordered derivatives.

In this study, we extend the application of PINNs to address inverse problems
involving time-fractional derivatives, specifically targeting two problems: 1)
anomalous diffusion and 2) fractional viscoelastic constitutive equation.
Leveraging both numerically generated datasets and experimental data, we
calibrate the concentration-dependent generalized diffusion coefficient and
parameters for the fractional Maxwell model. We devise a tailored residual
loss function that scales with the standard deviation of observed data.

We rigorously test our framework’s efficacy in handling anomalous diffusion.
Even after introducing 25% Gaussian noise to the concentration dataset, our
framework demonstrates remarkable robustness. Notably, the relative error in
predicting the generalized diffusion coefficient and the order of the
fractional derivative is less than 10% for all cases, underscoring the
resilience and accuracy of our approach. In another test case, we predict
relaxation moduli for three pig tissue samples, consistently achieving
relative errors below 10%. Furthermore, our framework exhibits promise in
modeling anomalous diffusion and non-linear fractional viscoelasticity.

##  1 Introduction

Fractional calculus broadens the traditional concepts of differentiation and
integration to include non-integer orders. Its roots trace back to 1695, when
Leibniz posed the question of defining the $\frac{1}{2}$th order derivative in
a letter to L’Hospital. Since then, numerous mathematicians have contributed
to its evolution, shaping its theoretical foundations. While fractional
calculus is often perceived as distinct from integer-order calculus, it
essentially serves as an extension, offering a broader and more natural
framework. Here, integer-order derivatives emerge as special cases within this
framework [1].

The practical relevance of fractional calculus has significantly increased
since the 1970s, when it gained attention from applied scientists and
engineers. This surge in interest stemmed from the recognition that fractional
calculus provides enhanced models for many intricate phenomena across various
domains. These include electromagnetism [2], rheology [3, 4], fluid mechanics
[5], and biology [6]. Fractional calculus excels in capturing the memory and
non-local effects inherent in these systems, aspects often inadequately
described by integer-order derivatives. Fractional models frequently entail
numerous parameters, such as fractional orders and coefficients, which are not
directly measurable or specifiable [7]. Addressing inverse problems involving
fractional-order derivatives presents notable challenges concerning stability,
uniqueness, and computational cost.

Physics-Informed Neural Networks (PINNs) [8] have emerged as powerful
computational frameworks, especially for inverse problems. PINNs have been
used to solve a vast range of problems in fluid mechanics [9, 10], rheology
[11, 12], and mass transport [13, 14]. Despite their broad applicability, most
implementations of PINNs primarily address problems involving integer-order
derivatives, typically computed via automatic differentiation. However, a
growing interest exists in extending PINNs to handle fractional order
derivatives, which offer a more nuanced understanding of complex phenomena.
One promising approach involves utilizing the finite difference method to
compute fractional order derivatives, thereby paving the way for the
development of fractional-PINNs [15]. This novel methodology has already
yielded promising results, facilitating the discovery of fractional order
models for turbulent flow [16] and epidemiological dynamics [17]. However, it
is worth noting that the L1 scheme [18], commonly employed to calculate time-
fractional derivatives in these contexts, is tailored exclusively for sub-
diffusion with derivatives within the range $(0,1)$, limiting its
applicability to capturing only integer-ordered time derivatives.

This study focuses on time-fractional diffusion equations characterizing
anomalous diffusion and the fractional Maxwell model governing viscoelastic
behavior. To tackle these complex phenomena, we employ a finite-difference
method [19] tailored for fractional orders within the range $[0,1]$. This
approach facilitates the calculation of fractional-order derivatives, of which
integer-order derivatives are special cases. An overview of fractional
derivatives, anomalous diffusion, and fractional viscoelasticity is provided
in Sections 2.1, 2.2, and 2.3, respectively. Subsequently, in section 2.4, we
elucidate the methodology devised for solving inverse problems associated with
the fractional diffusion equation and the fractional Maxwell model. Notably,
our approach ensures that all losses are weighted according to the inherent
scales present in the data. Lastly, the outcomes of our investigations into
inverse problems are meticulously presented and analyzed in Section 3. Through
this structured approach, we aim to offer insights into the intricate
interplay between fractional calculus and real-world phenomena, shedding light
on their underlying mechanisms and behaviors.

##  2 Problem setup and methodology

###  2.1 Fractional derivatives

While offering profound insights into complex phenomena, fractional calculus
introduces notable computational and analytical challenges. The non-local
characteristics of fractional derivatives render them computationally
intensive and more challenging to implement than their integer-order
counterparts. Moreover, navigating the realm of fractional calculus demands a
firm grasp of intricate mathematical concepts, including fractional order
operators, differential equations, calculus of variations, and fractional
processes. Thus, a robust mathematical foundation coupled with a nuanced
understanding of its applications is imperative. The Caputo fractional
derivative of order $\alpha$, for $0<\alpha<1$ for a function $f$ and time
$t$, is defined as [1]

$\frac{\partial^{\alpha}f(t)}{\partial
t^{\alpha}}=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}\frac{\partial
f(s)}{\partial s}\frac{ds}{(t-s)^{\alpha}},$ (1)

where $\Gamma$ is the gamma function. For $\alpha=0$, Eq. (1) corresponds to
the classical Helmholtz elliptic equation, and for the limit
$\alpha\rightarrow 1$, the integer order time derivative can be obtained. For
time-fractional differential equations, it is possible to numerically
approximate the value of the fractional derivative. For $0\leq\alpha\leq 1$,
the finite difference approximation of the Caputo fractional derivative
($L_{t}^{\alpha}$) of order $\alpha$ can be evaluated as[19]

$L_{t}^{\alpha}f\left(t_{k+1}\right):=\frac{1}{\Gamma(2-\alpha)}\sum_{j=0}^{k}b_{j}\frac{f(t_{k+1-j})-f(t_{k-j})}{\Delta
t^{\alpha}},$ (2)

where k is the time step number and $\Delta t$ is the time step size. The
coefficients $b_{j}$ are defined as

$b_{j}:=(j+1)^{1-\alpha}-j^{1-\alpha}.$ (3)

We rewrite this equation to get

$L_{t}^{\alpha}f\left(t_{k+1}\right):=\frac{1}{\Delta
t^{\alpha}\Gamma(2-\alpha)}(b_{0}f(t_{k+1})-b_{0}f(t_{k})+\sum_{j=0}^{k-1}b_{j+1}f(t_{k-j})-\sum_{j=1}^{k}b_{j}f(t_{k-j})).$
(4)

###  2.2 Anomalous diffusion

Diffusion, driven by the erratic thermal motion of molecules, underpins
numerous natural phenomena across physical, chemical, and biological realms.
In classical Fickian diffusion, transport dynamics align neatly with Gaussian
statistics, as evidenced by the linear time evolution of mean squared
displacement ($<x^{2}(t)>$). This linear relationship, captured by the
equation [20]:

$<x^{2}(t)>\sim D_{1}t,$ (5)

where $D_{1}$ is the diffusion coefficient illustrates that the average
distance traversed by diffusing particles grows proportionally with the square
root of time.

However, in intricate systems like porous media, amorphous semiconductors, and
polymeric structures, transport behavior often diverges from this simplistic
pattern, heralding what is known as anomalous diffusion. Anomalous diffusion
is typified by a non-linear augmentation of mean squared displacement over
time, characterized by the power-law expression:

$<x^{2}(t)>\sim\tilde{D}t^{\alpha}.$ (6)

Here, $\alpha$ denotes the anomalous diffusion exponent, and $\tilde{D}$
denotes the generalized diffusion coefficient. The exponent $\alpha$ measures
the degree of deviation from normal diffusion. When $\alpha=1$, the diffusion
is normal and follows Fick’s law. When $\alpha>1$, the diffusion is faster
than Fickian diffusion and is called super-diffusion. When $\alpha<1$, the
diffusion is slower than Fickian diffusion and is called sub-diffusion.
Typical diffusion processes can be described by the diffusion equation with
integer order derivatives

$\frac{\partial c}{\partial t}=\nabla\cdot(D_{1}\nabla c),$ (7)

Here, c denotes the concentration of diffusing particles and $\nabla$ denotes
the gradient operator. Processes with anomalous diffusion are described by
fractional order derivatives. For example, sub-diffusion can be modeled by
using a time-fractional derivative of order $\alpha$, where $0<\alpha<1$, as
follows

$\frac{\partial^{\alpha}c}{\partial t^{\alpha}}=\nabla\cdot(\tilde{D}\nabla
c).$ (8)

where $\tilde{D}$ denotes the sub-diffusion coefficient. This equation implies
that the diffusing particles have a long-tailed waiting time distribution,
meaning they tend to stay longer in some regions before moving to another. The
jump length distribution, however, remains finite. For a concentration-
dependent diffusion coefficient, we can write [21]

$\frac{\partial^{\alpha}c}{\partial
t^{\alpha}}=\tilde{D}(c_{xx}+c_{yy})+\tilde{D}_{c}(c_{x}^{2}+c_{y}^{2}),$ (9)

where the subscripts $x$, $y$ and $c$ denote the derivatives with respect to
the $x$-coordinate, $y$-coordinate and the concentration, respectively.

###  2.3 Fractional viscoelasticity

The theory of linear viscoelasticity offers a robust mathematical framework
for understanding the interplay between stress ($\tau(t)$), strain
($\epsilon(t)$), and time ($t$). A common examination method, known as a
relaxation test, involves subjecting a material to a constant strain and
analyzing its stress response. In such tests, the resulting stress ($\tau(t)$)
can be expressed as:

$\tau(t)=G(t)\epsilon^{0},$ (10)

where $G(t)$ is defined as the relaxation modulus and $\epsilon^{0}$ is the
strain applied at the initial time. The traditional linear visoelastic models
are built using springs and dashpots. An ideal spring, obeying Hooke’s law,
relates stress and strain as:

$\tau(t)=E\epsilon(t),$ (11)

while a dashpot, adhering to Newton’s model, equates stress to the derivative
of strain with respect to time:

$\tau(t)=\eta\frac{d\epsilon(t)}{dt},$ (12)

where E is the elastic modulus and $\eta$ is the material viscosity. These
elements can be interconnected to form more complex models, with the
relaxation modulus expressions derived accordingly. Connecting a spring and a
dashpot in series yields the Maxwell model, while parallel connection results
in the Kelvin-Voigt model. Introducing additional elements facilitates the
capture of multiple timescales in material behavior. However, the
incorporation of numerous elements poses challenges. Interpretation becomes
more intricate, and the model becomes increasingly susceptible to overfitting
as the number of elements rises. An alternative approach to capturing complex
material behavior with fewer parameters involves the utilization of a
viscoelastic element known as the springpot. In the case of a springpot, the
stress-strain relationship is defined as:

$\tau(t)=K_{\alpha}\frac{d^{\alpha}\epsilon(t)}{dt^{\alpha}},$ (13)

where $K_{\alpha}$ is the ‘firmness’ of the material and $\alpha$ is the order
of the fractional derivative. The springpot’s behavior varies with the value
of the fractional derivative order. The stress-strain relationship for a
fractional Maxwell model is given by [3]

$\tau(t)+\eta\frac{d^{\nu}\tau(t)}{dt^{\nu}}=\kappa\frac{d^{\nu}\epsilon(t)}{dt^{\nu}},$
(14)

with the relaxation modulus being defined as

$G(t)=\frac{\kappa}{\eta}E_{\nu}\left(-\frac{t^{\nu}}{\eta}\right),$ (15)

where $E_{\nu}$ is Mittag-Leffler function of order $\nu$. The Mittag-Leffler
function is defined as

$E_{\nu}(z)=\sum_{n=0}^{\infty}\frac{z^{n}}{\Gamma(\nu n+1)},$ (16)

where $\Gamma$ is the gamma function. The fractional Maxwell model has the
capability of capturing long memory effects and can capture complex material
behavior with fewer parameters.

###  2.4 Physics-Informed Neural Networks

Figure 1: The problem set up for (a) anomalous diffusion and (b) the
fractional Maxwell model, as described in section 2.4. Here, I denotes the
identity operator, and the differential operators $\partial x$, $\partial y$
and $\partial c$ are computed using automatic differentiation.

While automatic differentiation enables the computation of integer-order
derivatives, direct calculation of fractional-order derivatives presents a
challenge. To address this, in computing the residual of a time-fractional
differential equation, we employ automatic differentiation for spatial order
derivatives and finite differences for time derivatives [15], utilizing
equation (4).

In our pursuit to infer the anomalous diffusion coefficient as a function of
concentration and fractional order derivative, we approximate the functions
$(t,x,y)\mapsto(c^{pu})$ and $(c^{pu})\mapsto(D)$ using two neural networks
parameterized by $\theta$ and $\phi$, respectively. Here, the superscripts
‘$pu$’ denote the physics uninformed nature of these networks.

To quantify the loss for regression over the concentration field, we define
the mean squared loss as follows:

$L_{data-
AD}(\theta)=\mathbb{E}_{(t,x,y,c)}[\frac{|{c^{pu}}(t,x,y;\theta)-{c}|^{2}}{{\sigma_{c}}^{2}}],$
(17)

where $\mathbb{E}$ denotes the expectation approximated by the population
mean, ${\sigma_{c}}$ is the standard deviation of the reference concentration
field $c$ and the subscript $AD$ is used to denote anomalous diffusion. We now
define

$h(c,t^{k})=b_{0}c(t^{k})-\sum_{j=0}^{k-1}b_{j+1}c(t^{k-j})+\sum_{j=1}^{k}b_{j}c(t^{k-j}),$
(18)

where $t^{k}$ is the $k^{th}$ time step. Further, let

$g(c,D)={D}(c_{xx}+c_{yy})+{D}_{c}((c_{x})^{2}+(c_{y})^{2}),$ (19)

as defined in Eq. (9). To construct the physics-informed network for anomalous
diffusion, we combine the equations (4) and (9) to get

$c^{pi}(t^{k+1},\boldsymbol{x};\Delta
t,\theta,\phi,\alpha)=\Gamma(2-\alpha)\Delta
t^{\alpha}g(c^{pu}(t^{k},\boldsymbol{x};\theta),\tilde{D}(c^{pu};\phi))+h(c^{pu}(t^{k},\boldsymbol{x};\theta),t^{k}),$
(20)

where $\boldsymbol{x}=(x,y)$ are the spatial coordinates and the superscript
‘$pi$’ is used to denote ‘physics-informed’. Since the physics-uninfomed and
physics-informed concentration fields evaluated at the same spatio-temporal
point need to be consistent, this allows us to define a consistency loss [22]
as

$L_{consistency-AD}(\theta,\phi;\Delta
t,\alpha)=\mathbb{E}_{(t,\boldsymbol{x})}[\frac{|{c}^{pi}(t,\boldsymbol{x};\Delta
t,\theta,\phi,\alpha)-{c}^{pu}(t,\boldsymbol{x};\theta)|^{2}}{{\sigma_{c}}^{2}}].$
(21)

The parameters $\theta$ and $\phi$ can now be optimized by minimizing the
following mean squared error function

$L_{MSE-AD}(\theta,\phi;\alpha)=L_{data-AD}(\theta)+L_{consistency-
AD}(\theta,\phi;\alpha).$ (22)

Similarly, to learn the parameters of the fractional Maxwell equation (Eq.
(14)), we utilize data on stress ($\tau$) and strain ($\epsilon$) as functions
of time ($t$). We approximate the functions $(t)\mapsto(\sigma^{pu})$ and
$(t)\mapsto(\epsilon^{pu})$ using two fully-connected neural networks
parameterized by $\beta$ and $\zeta$, respectively. The mean squared loss for
regression over the stress and strain fields is defined as:

$L_{data-
FM}(\beta,\zeta)=\mathbb{E}_{(t,\tau)}[\frac{|{\tau^{pu}}(t;\beta)-{\tau}|^{2}}{{\sigma_{\tau}}^{2}}]+\mathbb{E}_{(t,\epsilon)}[\frac{|{\epsilon^{pu}}(t;\zeta)-{\epsilon}|^{2}}{{\sigma_{\epsilon}}^{2}}],$
(23)

where ${\sigma_{\tau}}$ and ${\sigma_{\epsilon}}$ are the standard deviations
of the reference stress field $\tau$ and the reference strain field,
respectively. To get the physics-informed stress, we combine equations (4) and
(14) to derive

$\tau^{pi}(t^{k+1};\Delta t,\beta,\zeta,\nu)=\kappa
L_{t}^{\nu}(\epsilon^{pu}(t^{k+1};\Delta t,\zeta))-\eta
L_{t}^{\nu}(\tau^{pu}(t^{k+1};\Delta t,\beta)),$ (24)

which allows us to define

$L_{consistency-FM}(\beta,\zeta;\Delta
t,\nu)=\mathbb{E}_{(t)}[\frac{|{\tau}^{pi}(t;\Delta
t,\beta,\zeta,\nu)-{\tau}^{pu}(t;\beta)|^{2}}{{\sigma_{\tau}}^{2}}],$ (25)

where the subscript $FM$ is used to denote fractional Maxwell. The parameters
$\beta$ and $\zeta$ are then optimized by minimizing the following mean
squared error function

$L_{MSE-FM}(\beta,\zeta)=L_{data-FM}(\beta,\zeta)+L_{consistency-
FM}(\beta,\zeta;\nu).$ (26)

##  3 Results and discussion

###  3.1 Numerical dataset - anomalous diffusion

We constructed a synthetic dataset to serve as a benchmark and evaluate the
efficacy of our framework. Our neural network architecture comprised a fully
connected network with 4 layers and 16 neurons per layer to approximate the
function mapping $(t,x,y)\rightarrow(c^{pu})$, and another fully connected
network with 2 layers and 4 neurons per layer to approximate the subsequent
mapping $(c^{pu})\rightarrow(D)$. Throughout this study, we employed the swish
activation function for all neural networks.

In our training process, we adopted a cosine learning rate schedule [23].
Specifically, we set the maximum learning rate ($\eta_{max}$) to
$2.5\text{e-}03$ and the minimum learning rate ($\eta_{min}$) to
$2.5\text{e-}06$, enabling us to compute the learning rate ($\eta$)
dynamically using the following formula:

$A=A_{min}+0.5(A_{max}-A_{min})\left(1+\cos\left(\frac{T_{cur}}{T_{max}}\pi\right)\right),$

where $A_{min}$ and $A_{max}$ are the minimum and maximum values of the
learning rate, respectively, $T_{cur}$ denotes the current time step and
$T_{max}$ represents the total number of time steps. The neural networks
underwent training for 100,000 iterations and the optimized parameters were
obtained by minimizing the loss in Eq. (22).

To assess the robustness of our framework, we introduced Gaussian noise to the
concentration field. The resulting predicted generalized diffusion
coefficients are visualized in Fig. 2, along with the relative errors in the
predicted generalized diffusion coefficient and fractional derivative order.
Our solver effectively captured both the generalized diffusion coefficient and
the fractional derivative order across all tested scenarios.

Figure 2: (a) The reference generalized diffusion coefficient compared to
different levels of added noise. The relative errors in the (b) generalized
diffusion coefficient and (c) fractional order for different levels of added
noise.

###  3.2 Experimental dataset - fractional Maxwell model

To evaluate our framework’s performance on experimental data to learn the
parameters of the fractional Maxwell model, we leverage the nonlinear stress-
strain and corresponding relaxation modulus data obtained from shear rheology
experiments conducted on minipig tissues by Mitra et al. [24]. We plot a
sample of the experimentally obtained stress and strain curves in Fig. 3. Our
neural network architecture consists of a fully connected network with 2
layers and 20 neurons per layer, tasked with approximating the function
mapping $(t)\rightarrow(\sigma^{pu})$. Additionally, another fully connected
network with the same architecture is employed to approximate the subsequent
mapping $(t)\rightarrow(\epsilon^{pu})$.

We iteratively learn the parameters of equation (14) by minimizing the loss
function defined in equation (26). Subsequently, we utilize these learned
parameters to predict the relaxation modulus as defined in equation (15). Our
study encompasses three distinct tissue samples: neck, belly, and breast
tissue. The outcomes of our analysis are visualized in Fig. 4. Additionally,
we provide detailed numerical values for the relative errors in Table 1.
Notably, all relative errors for the predicted relaxation moduli concerning
the reported values are found to be less than 10%.

Figure 3: Plot of the experimentally obtained (a) stress and (b) strain values for a sample reported in [24].  Figure 4: The predicted relaxation modulus and the experimental data corresponding to three samples of minipig skin tissue.  Table 1: The relative errors for the relaxation moduli (G(t)) for the three pig tissue samples. Samples | Sample 1 | Sample 2 | Sample 3
---|---|---|---
Relative error | $9.14\times 10^{-2}$ | $3.06\times 10^{-2}$ | $8.92\times 10^{-2}$

##  4 Conclusions and future work

Time-fractional differential equations find widespread applications across
various fields, albeit often posing computational challenges and instability
issues, especially in the context of inverse problems. Recognizing this, there
is growing interest in leveraging Physics-Informed neural network-based
frameworks to address inverse problems entailing time-fractional derivatives.
In such frameworks, the fractional time derivative can be efficiently computed
using finite differences, while automatic differentiation handles other
derivatives.

In this study, we tackle two inverse problems associated with 1) anomalous
diffusion and 2) fractional viscoelasticity employing the aforementioned
approach. Our methodology involves defining the residual loss in a manner that
facilitates scaling the loss terms with the standard deviation of the observed
data. We utilize numerically generated datasets and experimental data for
learning the fractional coefficient and the concentration-dependent
generalized diffusion coefficient, as well as for calibrating parameters for
the fractional Maxwell model using stress and strain data over time.

To evaluate the performance of our framework in handling anomalous diffusion,
we conducted rigorous testing by comparing our predictions for a
concentration-dependent generalized diffusion coefficient and the fractional
order of the derivative against synthetically generated reference values.
Remarkably, even with the addition of a 25% Gaussian noise to the
concentration dataset, our framework demonstrated robust performance.
Specifically, we observed that the relative error in predicting the
generalized diffusion coefficient and the order of the fractional derivative
was less than 10%. This outcome underscores the resilience and accuracy of our
framework, even in the presence of significant noise, reaffirming its
reliability in practical applications.

To validate our findings, we predict relaxation moduli for three distinct
samples of pig tissue and compare them against reported values in the
literature. The relative errors are consistently below 10%, highlighting the
efficacy of the fractional model which requires fewer parameters. Furthermore,
our approach holds promise for extending to model non-linear fractional
viscoelasticity and can readily incorporate experimental data to validate
anomalous diffusion predictions. Additionally, the solver’s applicability can
be extended to three-dimensional scenarios, offering a broader scope for
future exploration into equations involving time-fractional derivatives.

## References

  * [1] Igor Podlubny.  Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications.  Elsevier, 1998.
  * [2] Nader Engheta.  On fractional calculus and fractional multipoles in electromagnetism.  IEEE Transactions on Antennas and Propagation, 44(4):554–566, 4 1996.
  * [3] Francesco Mainardi.  Fractional Viscoelastic Models.  Fractional Calculus and Waves in Linear Viscoelasticity, pages 57–76, 2010.
  * [4] A. Bonfanti, J. L. Kaplan, G. Charras, and A. Kabla.  Fractional viscoelastic models for power-law materials.  Soft Matter, 16(26):6002–6020, 2020.
  * [5] Zaid Odibat and Shaher Momani.  The variational iteration method: An efficient scheme for handling fractional partial differential equations in fluid mechanics.  Computers and Mathematics with Applications, 58(11-12):2199–2208, 2009.
  * [6] Anna Stankiewicz.  Fractional Maxwell model of viscoelastic biological materials.  BIO Web of Conferences, 10:02032, 2018.
  * [7] Bangti Jin and William Rundell.  A tutorial on inverse problems for anomalous diffusion processes.  Inverse Problems, 31(3), 2015.
  * [8] Maziar Raissi, Paris Perdikaris, and George Em 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.
  * [9] Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Karniadakis.  Physics-informed neural networks (PINNs) for fluid mechanics: a review.  Acta Mechanica Sinica/Lixue Xuebao, 37(12):1727–1738, 2021.
  * [10] Xiaowei Jin, Shengze Cai, Hui Li, and George Em Karniadakis.  NSFnets (Navier-Stokes Flow nets): Physics-informed neural networks for the incompressible Navier-Stokes equations.  (Hui Li).
  * [11] Mohammadamin Mahmoudabadbozchelou, George Em Karniadakis, and Safa Jamali.  nn-PINNs: Non-Newtonian physics-informed neural networks for complex fluid modeling.  Soft Matter, 18(1):172–185, 2022.
  * [12] Sukirt Thakur, Maziar Raissi, and Arezoo M. Ardekani.  ViscoelasticNet: A physics informed neural network framework for stress discovery and model selection.  pages 1–21, 2022.
  * [13] Jiangyu Wang, Xingjie Peng, Zhang Chen, Bingyan Zhou, Yajin Zhou, and Nan Zhou.  Surrogate modeling for neutron diffusion problems based on conservative physics-informed neural networks with boundary conditions enforcement.  Annals of Nuclear Energy, 176:109234, 2022.
  * [14] Sukirt Thakur, Ehsan Esmaili, Sarah Libring, Luis Solorio, and Arezoo M. Ardekani.  Inverse resolution of spatially varying diffusion coefficient using Physics-Informed neural networks.  pages 1–15, 2024.
  * [15] Guofei Pang, Lu Lu, and George Em Karniadakis.  fPINNs: Fractional Physics-Informed Neural Networks.  SIAM Journal on Scientific Computing, 41(4):A2603–A2626, 1 2019.
  * [16] Pavan Pranjivan Mehta, Guofei Pang, Fangying Song, and George Em Karniadakis.  Discovering a universal variable-order fractional model for turbulent couette flow using a physics-informed neural network.  Fractional Calculus and Applied Analysis, 22(6):1675–1688, 2019.
  * [17] Ehsan Kharazmi, Min Cai, Xiaoning Zheng, Zhen Zhang, Guang Lin, and George Em Karniadakis.  Identifiability and predictability of integer- and fractional-order epidemiological models using physics-informed neural networks.  Nature Computational Science, 1(11):744–753, 2021.
  * [18] Bangti Jin, Raytcho Lazarov, and Zhi Zhou.  An analysis of the L1 scheme for the subdiffusion equation with nonsmooth data.  IMA Journal of Numerical Analysis, 36(1):197–221, 2014.
  * [19] Yumin Lin and Chuanju Xu.  Finite difference/spectral approximations for the time-fractional diffusion equation.  Journal of Computational Physics, 225(2):1533–1552, 8 2007.
  * [20] Ralf Metzler and Joseph Klafter.  The random walk’s guide to anomalous diffusion: A fractional dynamics approach.  Physics Report, 339(1):1–77, 2000.
  * [21] Adib Ahmadzadegan, Jiacheng Zhang, Arezoo M. Ardekani, and Pavlos P. Vlachos.  Spatiotemporal measurement of concentration-dependent diffusion coefficient.  Physics of Fluids, 34(5), 2022.
  * [22] Sukirt Thakur, Maziar Raissi, Harsa Mitra, and Arezoo Ardekani.  Temporal Consistency Loss for Physics-Informed Neural Networks.  pages 1–14, 2023.
  * [23] Ilya Loshchilov and Frank Hutter.  SGDR: Stochastic gradient descent with warm restarts.  5th International Conference on Learning Representations, ICLR 2017 - Conference Track Proceedings, pages 1–16, 2017.
  * [24] Harsa Mitra, Evelyn Nonamaker, Ria D Cordera, Luis Solorio, and Arezoo M. Ardekani.  Rheological and Lipid Characterization of Minipig and Human Skin Tissue : A Comparative Study Across Different Locations and Depths .  2024\.


# Serial Interference Cancellation for Improving uplink in LoRa-like Networks

Angesom Ataklity TESFAY1, Eric Pierre SIMON1, Guillaume FERRÉ2 and Laurent
CLAVIER1,3
1Univ. Lille, CNRS, UMR 8520 - IEMN, F-59000, Lille, France (e-mail:
<EMAIL_ADDRESS>
2Univ. Bordeaux, CNRS, Bordeaux INP, IMS, UMR 5218, F-33400, Talence, France
(e-mail<EMAIL_ADDRESS>
3IMT Lille Douai, France (e-mail<EMAIL_ADDRESS>
###### Abstract

In this paper, we present a new receiver design, which significantly improves
performance in the Internet of Things networks such as LoRa, i.e., having a
chirp spread spectrum modulation. The proposed receiver is able to demodulate
multiple users simultaneously transmitted over the same frequency channel with
the same spreading factor. From a non-orthogonal multiple access point of
view, it is based on the power domain and uses serial interference
cancellation. Simulation results show that the receiver allows a significant
increase in the number of connected devices in the network.

###### Index Terms:

LoRa, SIC, IoT, CSS

##  I Introduction

Internet of Things (IoT) offers interconnection between a large number of
physical objects used for different applications such as environmental
monitoring, smart farming and metering. The number of devices in the IoT
network is expected to increase exponentially [1, 2] in the coming years. Over
recent years, Low Power Wide Area Networks (LPWAN) technologies, like LoRa and
SigFox in unlicensed bands, are proposed to face this growth. LoRa employs a
Chirp-Spread Spectrum (CSS) Modulation technique to encode information [5, 6].
Several Spreading Factors (SF) and bandwidth settings enable quasi-orthogonal
transmissions for LoRa end-devices. However, when two or more end-devices
transmit simultaneously on the same band and with the same SF, a collision
occurs at the receiver. This collision generally results in the loss of all
colliding packets, but sometimes the first arrived signal might be fully
decoded due to capture effect.

The authors in [7] examines the interfering signal and analytically derives
the performance of LoRa under same SF interference. Information of interferers
is lost in the process. Recent work [8] proposes an algorithm to decode the
correct symbols of synchronized or slightly desynchronized signals using their
timing information. However, ensuring synchronization is complex in the up-
link of long range networks.

Recently, the authors in [9] proposed a receiver which can decode two signals
received at the same time, with the same spreading factor using serial
interference cancellation (SIC). Asynchronicity is processed by estimating the
time shift between the two received signals. However, the authors presumed
that the signal with highest power was received first and this is not
necessarily the case. The extension of this work, which is able to process
multiple simultaneously received signals with the same SF is presented in
[10]. However, the power ratio between two successive received signals is
assumed to be constant.

The work in [11], closely related to ours, proposes SIC to decode information
from all users. A synchronized reception is considered, which is not realistic
in the uplink.

In this work, we propose a SIC technique to enable a receiver to decode
simultaneously multiple signals. We extend the work in [11] to a more general
case where transmitters are asynchronous. Our main contributions are to
develop a complete receiver structure, including the detection of packets,
channel estimation, detection of symbols and interference cancellation. We
analyse the performance in the presence of multiple interfering nodes and show
that this approach can significantly improve the scalability of the network.
This work is divided into five sections. Section II describes some key
features of LoRa and introduces the system model. Section III presents the SIC
algorithm while section IV discusses the simulation results. Conclusions are
drawn in section V.

$\displaystyle
s_{p,1}^{(j,i)}(t)=\begin{cases}&\exp\bigg{(}2\jmath\pi\Big{(}\frac{B}{2T_{s}}t^{2}+\big{(}\frac{m_{p,1}^{(i)}-B\Delta^{(j,i)}}{T_{s}}\big{)}t+\phi_{p,1}^{(j,i)}\Big{)}\bigg{)},\\\
&\hskip 142.26378ptt\in
A_{p,1}=\big{[}\frac{-T_{s}}{2}+\Delta^{(j,i)},\frac{T_{s}}{2}+\Delta^{(j,i)}-\tau_{p,1}^{(i)}\big{[},\\\
&\exp\bigg{(}2\jmath\pi\Big{(}\frac{B}{2T_{s}}t^{2}+\big{(}\frac{m_{p,1}^{(i)}-B\Delta^{(j,i)}}{T_{s}}\big{)}t+\phi_{p,1}^{(j,i)}-B(t-\Delta^{(j,i)})\Big{)}\bigg{)},\\\
&\hskip 142.26378ptt\in
B_{p,1}=\big{[}\frac{T_{s}}{2}+\Delta^{(j,i)}-\tau_{p,1}^{(i)},\frac{T_{s}}{2}\big{[}.\end{cases}$
(4)

$\displaystyle
s_{p,2}^{(j,i)}(t)=\begin{cases}&\exp\bigg{(}{2\jmath\pi\Big{(}\frac{B}{2T_{s}}t^{2}+\big{(}\frac{m^{(i)}_{p,2}+B(T_{s}-\Delta^{(j,i)})}{T_{s}}\big{)}t+\phi_{p,2}^{(j,i)}\Big{)}\bigg{)}},\\\
&\hskip 142.26378ptt\in
B_{p,2}=\big{[}\frac{-T_{s}}{2},\frac{-T_{s}}{2}+\Delta^{(j,i)}-\tau_{p,2}^{(i)}\big{[},\\\
&\exp\bigg{(}2\jmath\pi\Big{(}\frac{B}{2T_{s}}t^{2}+\big{(}\frac{m_{p,2}^{(i)}+B(T_{s}-\Delta^{(j,i)})}{T_{s}}\big{)}t+\phi_{p,2}^{(j,i)}-B(t+T_{s}-\Delta^{(j,i)})\Big{)}\bigg{)},\\\
&\hskip 142.26378ptt\in
A_{p,2}=\big{[}\frac{-T_{s}}{2}+\Delta^{(j,i)}-\tau^{(i)}_{p,2},\frac{-T_{s}}{2}+\Delta^{(j,i)}\big{[}.\end{cases}$
(5)

$\text{ where
}\;\;\;\;\;\phi_{p,1}^{(j,i)}=\frac{B(\Delta^{(j,i)})^{2}-2m_{p,1}^{(i)}\Delta^{(j,i)}}{2T_{s}},\;\;\text{and}\;\;\phi_{p,2}^{(j,i)}=\frac{B(T_{s}-\Delta^{(j,i)})^{2}+2m_{p,2}^{(i)}(T_{s}-\Delta^{(j,i)})}{2T_{s}}.$

##  II System Model

The LoRa signal model is an $M$-ary digital modulation based on frequency
shift chirp spread spectrum. The $M$ possible waveforms are chirps modulated
over the frequency band $B$ with $M$ different initial frequencies. A
modulating symbol carries $\textrm{SF}=\log_{2}(M)$ bits, with SF ranging from
7 to 12111$\textrm{SF}=6$ also exists but the modulation scheme is modified..
The symbol duration is $T_{s}=MT$, where $T=1/B$ is the sampling period. Each
modulated chirp is based on a raw chirp, whose instantaneous frequency is
given by $f(t)=\frac{B}{T_{s}}t$. Recall that the instantaneous phase
($\vartheta(t)$: $f(t)=\frac{1}{2\pi}\frac{d\vartheta(t)}{dt}$) is the
integral of $f(t)$, which yields the following baseband expression for the raw
chirp:

$c(t)=\exp\Big{(}2\jmath\pi\frac{B}{T_{s}}\frac{t^{2}}{2}\Big{)},\;t\in\Big{[}-\frac{T_{s}}{2},\frac{T_{s}}{2}\Big{[}.$
(1)

The transmitted modulating symbol of the $i$th user at time
$pT_{s},\;p=0,\cdots,P-1$, with $P$ the number of symbols transmitted in a
packet, is denoted by $m^{(i)}_{p}\in\\{0,...,M-1\\}$. The corresponding
modulated chirp is obtained by left-shifting the raw chirp of
$\tau^{(i)}_{p}=m^{(i)}_{p}T$ in the time domain, yielding a shift of
$\frac{m^{(i)}_{p}}{T_{s}}$ in the frequency domain (see Fig. 1):

$\displaystyle
s_{p}^{(i)}(t)=\begin{cases}&\exp\big{(}2\jmath\pi\big{[}\frac{B}{2T_{s}}t^{2}+\frac{m_{p}^{(i)}}{T_{s}}t\big{]}\big{)},\\\
&\hskip
71.13188ptt\in\text{\big{[}$-\frac{T_{s}}{2},\frac{T_{s}}{2}-\tau^{(i)}_{p}$}\big{[},\\\
&\exp\big{(}2\jmath\pi\big{[}\frac{B}{2T_{s}}t^{2}+\big{(}\frac{m_{p}^{(i)}}{T_{s}}-B\big{)}t\big{]}\big{)},\\\
&\hskip 71.13188ptt\in\text{
\big{[}$\frac{T_{s}}{2}-\tau^{(i)}_{p},\frac{T_{s}}{2}$}\big{[}.\end{cases}$
(2)

Then, the complex envelope of the transmitted signal of user $i$,
$s^{(i)}(t)$, is given by:

$s^{(i)}(t)=\sum_{p=0}^{P-1}s_{p}^{(i)}(t-pT_{s}).$ (3)

The transmitted LoRa packet is finally obtained by adding the preamble
consisting of $N_{p}$ raw chirps.

Figure 1: Raw up-chirp and Coded chirp associated with $m^{(i)}_{p}$.

The rest of this section deals with the reception of a LoRa frame in the
presence of interfering users. Indeed, when two users simultaneously transmit
their data using the same spreading factor and frequency band, a collision can
possibly occur at the reception. Without loss of generality, only one
interfering user is considered. Complete expressions with $N_{u}$ interferers
are directly obtained by adding all interfering terms, which all have a
similar expression.

We describe now the collision between a $p$th symbol of a user of interest $j$
and two consecutive symbols of an interfering user $i$ as shown in Fig. 2.
Recall that $s^{(j)}_{p}(t)$ denotes the modulated symbol (or chirp) at time
$pT_{s}$ of user $j$ (see (2)). The two consecutive interfering symbols of
user $i$ are denoted by $s_{p,1}^{(j,i)}(t)$ and $s_{p,2}^{(j,i)}(t)$, the
index $j$ indicating that the synchronization is carried out on user $j$. The
shifts carrying the corresponding modulating symbols are given by
$\tau^{(i)}_{p,1}=m_{p,1}^{(i)}T$ and $\tau^{(i)}_{p,2}=m_{p,2}^{(i)}T$. The
delay between user $j$ and user $i$ (accounting for asynchronous
transmissions) is denoted by $\Delta^{(j,i)}$. To formulate the mathematical
expression for the interfering signals in Fig. 2, we consider two conditions
that define the positions of the shifts of the interfering symbols:

$\displaystyle\
\text{$C_{1}$:}\;\tau^{(i)}_{p,1}>\Delta^{(j,i)},\text{$C_{2}$:}\;\;\;\tau^{(i)}_{p,2}<\Delta^{(j,i)}.$
(6)

Depending on these conditions, there could be two to four contributions to the
interference.

The expression for the interfering signal is given in (4) and (5). The second
case in (4) (respectively the first in (5)) contributes only if $C_{1}$
(respectively $C_{2}$) in (6) is satisfied.

Figure 2: LoRa-like Interfering symbols.

To synchronize, we perform a correlation between the received signal and the
preamble. The highest peak indicates the strongest user, indexed by $j$, that
we will try to decode first, so we synchronize on it.

The received signal associated with symbol $p$ of user $j$ sampled at $t=nT$,
$n=-\frac{M}{2},\cdots,\frac{M}{2}-1$ is:

$\displaystyle y_{p}^{(j)}[n]=h^{(j)}$
$\displaystyle\sqrt{P_{t}}\;s_{p}^{(j)}[n]$
$\displaystyle+h^{(i)}\sqrt{P_{t}}\left(s_{p,1}^{(j,i)}[n]+s_{p,2}^{(j,i)}[n]\right)+w_{p}[n].$
(7)

where we have used the notation $s[n]=s(nT_{s})$ for any modulated signal
$s(t)$. $w_{p}[n]\sim\mathcal{CN}(0,\,\sigma_{n}^{2})$ is a complex Gaussian
noise, $h^{(j)}$ and $h^{(i)}$ are the channel attenuations for users $j$ and
$i$, $P_{t}$ is the transmitted power (we consider the same transmit power for
each user).

$\displaystyle r_{p}[n]$
$\displaystyle=h^{(j)}\sqrt{P_{t}}\;e^{2\jmath\pi\frac{m_{p}^{(j)}}{M}n}$
$\displaystyle+h^{(i)}\sqrt{P_{t}}\bigg{(}\sum_{l=1}^{2}e^{2\jmath\pi\phi_{p,l}^{(j,i)}}e^{2\jmath\pi\big{(}\frac{\varphi_{p,l}^{(j,i)}}{M}\big{)}n}\mathds{1}_{A_{p,l}}[n]+e^{2\jmath\pi(\phi_{p,l}^{(j,i)}+B\Delta^{(j,i)})}e^{2\jmath\pi\big{(}\frac{\varphi_{p,l}^{(j,i)}}{M}-(-1)^{l+1}\big{)}n}\mathds{1}_{B_{p,l}}[n]\bigg{)}+w_{p}[n],$
(8)

$R^{\prime}_{p}[k]=P^{(j)}\delta[k-m_{p}^{(j)}]+P^{(i)}\Bigg{(}\sum_{l=1}^{2}e^{\jmath\beta^{(j,i)}_{p,l}}\frac{\sin\big{(}\pi\big{(}\frac{k-\varphi^{(j,i)}_{p,l}}{M}\big{)}\Delta
A_{p,l}\big{)}}{\sin\big{(}\pi\big{(}\frac{k-\varphi^{(j,i)}_{p,l}}{M}\big{)}\big{)}}+e^{\jmath\gamma^{(j,i)}_{p,l}}\frac{\sin\big{(}\pi\big{(}\frac{k-\varphi^{(j,i)}_{p,l}}{M}-(-1)^{l}\big{)}\Delta
B_{p,l}}{\sin\big{(}\pi\big{(}\frac{k-\varphi^{(j,i)}_{p,l}}{M}-(-1)^{l}\big{)}\big{)}}\Bigg{)}+W^{\prime}_{p}[k].$
(9)

$\displaystyle\text{where}\;\beta^{(j,i)}_{p,l}=2\pi\phi_{p,l}^{(j,i)}-\pi\big{(}\frac{k-\varphi^{(j,i)}_{p,l}}{M}\big{)}(\Delta
A_{p,l}-1),\;\gamma^{(j,i)}_{p,l}=2\pi(\phi_{p,l}^{(j,i)}+B\Delta^{(j,i)})-\pi\Big{(}\frac{k-\varphi^{(j,i)}_{p,l}}{M}-(-1)^{l}\Big{)}(\Delta
B_{p,l}-1),$
$\displaystyle\varphi_{p,l}^{(j,i)}={m_{p,l}^{(i)}-B\Delta^{(j,i)}},\;\Delta
A_{p,l},\;\Delta B_{p,l}\;\text{are the length of the
intervals}\;A_{p,l}\;\text{and}\;B_{p,l}\;\text{for}\;l\in\\{1,2\\}.$

To demodulate, the samples of the received signal are multiplied by the
conjugate of the raw chirp: $r_{p}[n]=y_{p}^{(j)}[n]c^{*}[n]$, where
$c^{*}[n]=c^{*}(nT)$. The resulting signal is given in (II), with $A_{p,l}$
and $B_{p,l}$, $l\in\\{1,2\\}$, from (4) and (5).

Taking the absolute value of the Fast Fourier Transform (FFT) of (II), the
useful information creates a peak whose amplitude depends on the channel
attenuation. Eq. (II) also shows that user $i$ can create up to four peaks in
the Fourier domain (and this is the same for any interfering user). The height
of these peaks is determined not only by the channel attenuation but also by
the length of the overlaps, and this is good news for our SIC scheme. Fig. 3
illustrates this phenomenon by showing the results of the FFT processing for
two coded chirps with $m=4$ and $m=99$ for the user of interest. The
parameters are as follows: Noise level of $-114$ dBm, $\textrm{SF}=8$, $B=250$
kHz and one interfering user.

##  III Proposed SIC Receiver and Detection

We consider SIC as a pure receiver technique, without any power control or any
types of adjustment at the transmitter side. First, the strongest signal is
decoded. Its contribution can then be reconstructed, using an estimate of the
channel, and subtracted from the received signal. From the resulting residue,
the second strongest signal can then be extracted [12]. The performance of the
decoding will be degraded by two facts: if the difference in received powers
from two different users is too close (which means that their channel is
similar because we do not perform any kind of power control); and if other
peaks collide making the peak of interest not the strongest one and/or
significantly deviating from its expected amplitude. It is important to notice
in this case that the signal to interference plus noise power ratio, with an
interference power evaluated through the sum of all interfering signal powers,
is not a reliable metric. We can also note that the wide range of the
communications should allow a large dispersion in the received powers,
allowing our approach to perform well. The decoding of a user can be described
as follows. First, the strongest preamble is detected by calculating the
correlation between the received signal and the preamble and taking the
maximal value, denoted by $C_{\textrm{max}}$. This allows us to detect the
beginning of the packet of the strongest user and to estimate its channel
$\hat{h}^{(j)}$. Multiplication by the down chirp is then performed, see (II).
We apply an FFT to $r_{p}[n]$ and divide it by the estimated channel, which is
shown in $R^{\prime}_{p}[k]={R_{p}[k]}/{\hat{h}^{(j)}}$,
$P^{(j)}=\sqrt{P_{t}}\;h^{(j)}/\hat{h}^{(j)}$,
$P^{(i)}=\sqrt{P_{t}}\;h^{(i)}/\hat{h}^{(j)}$,
$W^{\prime}_{p}[k]={W_{p}[k]}/{\hat{h}^{(j)}}$, and $W_{p}[k]$ is the FFT of
the Gaussian noise. In (9), it is shown that the received power varies
according to the standard deviation $\sigma_{n}$ of the noise term
$W^{\prime}_{p}[k]$. (Expressions are still given for two users and are simply
extended to $N_{u}$ users.)

The next step is to search for the peaks of the desired user. Instead of the
usual maximum search, we propose to search for the peak having the closest
value to $P^{(j)}$ by minimizing the difference between the detected peaks in
$R^{\prime}_{p}[k]$ and $P^{(j)}$. This avoids mistakes when interferers
collide and create peaks with a larger amplitude than the desired one. We also
check that $|R^{\prime}_{p}[k]-P^{(j)}|<\varepsilon$ to detect collisions on
the desired peak, where $\varepsilon$ is tuned by grid-search, yielding an
optimal value of $3\sigma_{n}$. If the above condition is not verified we
simply choose the maximum peak. Once the symbols of the desired user is
decoded, its contribution is re-constructed and subtracted from the composite
received discrete signal $y_{p}[n]$. The process then starts again to detect
the next strongest user. This detection/suppression process is repeated until
no more preamble can be found or a packet is not correctly decoded. Algorithm
1 illustrates the proposed receiver. The input is only the combined received
signal.

Algorithm 1 SIC Receiver for LoRa-like Networks

0: Received samples $y_{p}[n]$, number of transmitted symbols $N$, detection
threshold $\varepsilon$

0: Decoded symbols $m^{(j)}_{p}$

1: $Detected\leftarrow 1$

2: while $Detected==1$ do

3: if $|C_{\textrm{max}}|<\varepsilon$ then

4: $Detected\leftarrow 0$

5: else

6: Estimate $\hat{h}^{(j)}$

7: $r_{p}[n]\leftarrow y_{p}^{(j)}[n]c^{*}[n]$

8: $R^{\prime}_{p}[k]\leftarrow|FFT(r_{p}[n])|$

9: for $p=1$ to $P$ do

10: $u\leftarrow\operatorname*{argmin}\limits_{k}|R^{\prime}_{p}[k]-P^{(j)}|$

11: if $|R^{\prime}_{p}[u]-P^{(j)}|^{2}<\varepsilon$ then

12: $m_{p}^{(j)}\leftarrow u$

13: else

14:
$m_{p}^{(j)}\leftarrow\operatorname*{argmax}\limits_{k}(|R^{\prime}_{p}[k]^{2}|)$

15: end if

16: end for

17: $\hat{s}_{p}^{(j)}[n]\leftarrow$ CSS Modulation of $m_{p}^{(j)}$

18: $y_{p}[n]\leftarrow y_{p}[n]-\hat{h}^{(j)}\;\hat{s}_{p}^{(j)}[n]$

19: end if

20: return $m_{p}^{(j)}$

21: end while

##  IV Simulation Result

###  IV-A Simulation set-up

In order to evaluate the performance of our proposed scheme, we rely on Monte
Carlo simulations of 1000 rounds. We consider a circle of radius $R$ with one
gateway at the center. Multiple users are uniformly distributed within the
circle. The distance from user $i$ to the gateway is denoted by $d^{(i)}$. We
consider a set of users transmitting their data through block fading channels
during the time interval of interest. All users employ CSS modulation with a
symbol duration $T_{s}$ (same SF). They transmit their signals with the same
power $P_{t}$, i.e., there is no power control. We consider asynchronous
transmission among the nodes, different devices operating autonomously. We
consider path loss and Rayleigh multi-path fading $\chi_{i}$. For path loss,
the signal amplitude decays with distance according to ${d^{(i)}}^{-\eta/2}$,
where $\eta$ is the path loss exponent. The channel attenuation (in amplitude)
is then $h^{(i)}={d^{(i)}}^{-\eta/2}\chi_{i}$. In the following we take a
maximum range $R=5$ km and we choose $\eta=3.5$. However, the channel
attenuation has to be such that the user can be connected to the network with
the chosen SF, in other words if their received power is greater than the
receiver sensitivity $S$, for instance $S=-124$ dBm for $\textrm{SF}=8$ and
$B=250$ kHz. Users that do not respect this condition are discarded and drawn
again. The noise level of a receiver at room temperature is
$-174+10\log_{10}(B)+\textrm{NF}=-114$ dBm, where NF is the receiver noise
figure and a classic $6$ dB is considered [13].

Figure 3:  FFT of coded chirp of the user of interest at the values of 4 and
99 in the presence of one interfering user using $\textrm{SF}=8$, $B=250$ KHz.

We consider a window of size $\textrm{W}=3T_{f}$, where $T_{f}$ is the frame
length. $N_{u}$ devices are transmitting in this window and asynchronicity is
ensured through a random variable $\Delta^{(i)}$, uniformly distributed over
$\left[0,2T_{f}\right]$. Indeed, $2T_{f}$ is the vulnerability period of an
ALOHA protocol and to be sure we receive all full packets, we need to end our
study frame at $3T_{f}$. A perfect time division would allow 3 users (3
packets) in this frame if collisions are not allowed.

No channel coding is performed and we are interested in the symbol error rate
(SER).

###  IV-B Performance for a given useful link.

For a given useful link characterized by its SNR
$\frac{|h|^{2}P}{\sigma_{n}^{2}}$, Fig. 4 shows the performance of this
specific receiver when there is no other interfering user ($N_{u}=1$) as well
as when $N_{u}=4$, $7$, and $10$ interfering users are present. The positions
and channels of these interfering users are randomly chosen at each round.

Figure 4: SER for $N_{u}=1$, $4$, $7$, and $10$, with $\textrm{SF}=8$, $B=250$
kHz, and different SNR values.

For an SNR of $-10$ dB, a typical value for $\textrm{SF}=8$ and $B=250$ kHz,
the performance is good even for 10 interfering users, we reach a SER of
$10^{-3}$ at SNR of $-10$ dB which is a target SER for LoRa under same SF
interference [7].

###  IV-C Mean performance.

We now evaluate the performance for any user, what ever their channel (so
their SNR) is. Fig. 5 shows the average symbol error rate of the proposed
receiver with $\textrm{SF}=8$ when the position of all users are randomly
drawn at each round and when different noise levels are considered. It can be
seen that instead of the three users that could be handled with a perfect time
division multiple access (TDMA), which is even more than what an ALOHA
protocol as the one used in LoRa can handle, a SIC approach allows to have
about 20 users with a SER under $10^{-2}$.

Figure 5: SER for different values of $N_{u}$ with $\textrm{SF}=8$, $B=250$
kHz, and several noise levels. Figure 6: SER of each users for $N_{u}=10$,
$\textrm{SF}=8$, $B=250$ kHz, Noise level of $-114$ dBm.

SIC receiver can be affected by error propagation, consequence of the serial
decoding, reconstruction and suppression process. The results is a higher SER
for users with lower power. Fig. 6 illustrates the SER of each users when
$N_{u}=10$. The users are ordered according to their received power, from the
highest (1) to the lowest (10), confirming the error propagation.

###  IV-D SF orthogonality.

It was shown in [14] that SFs have imperfect orthogonality. Fig. 7 presents
the average SER for 10 users, which transmit their signal using
$\textrm{SF}=8$. We consider an additional number $N_{u}$ of interfering users
that transmit their signal using a different SF ($9,10$ or $12$). As long as
these added

Figure 7: SER of Fixed $10$ users using $\textrm{SF}=8$ and additional
interfering users using different SF (9, 10, 12), $B=250$ kHz, Noise level of
$-114$ dBm.

interferers are less than four, their impact is negligible. However, the
spreading factors are not perfectly orthogonal so that the performance of the
10 users using $\textrm{SF}=8$ degrades when the number of interfering users
with a different SF gets larger. We also notice that the shorter SF (9)
impacts more than the longer ones.

##  V Conclusion

Scalability is a major concern in IoT networks and actual deployed solutions
cannot face it. However, in the uplink communication, due to the star topology
of LPWANs, the complexity can be added at the receiver to improve the system
performance, without modifying the low consumption transmission scheme. And
this allows us to propose a serial interference cancellation scheme for chirp
spread spectrum based physical layer. Assuming uncoordinated transmitting
devices, we use a SIC from the power domain NOMA idea to recover multiple
signals arriving simultaneously, on the same bandwidth and with the same SF.
The diversity in received powers is ensured by the random radio channels faced
by the devices. A SIC receiver allows then to recover the multiple signals.

In fact, though the simulation results we showed that the proposed solution
allows supporting 20 times more devices without any modification of the
transmission scheme. This is attained with a first simple implementation of
the SIC and further improvement can certainly be obtained with a judicious
choice of the PHY layer parameters.

## Acknowledgement

This work was supported by IRCICA, USR CNRS 3380, Lille, by the COST action
CA15104, IRACON.

## References

  * [1] N. Shahid and S Aneja, _”Internet of Things: Vision, application areas and research challenges,”_ Int. Conf. on I-SMAC, Palladam, 2017.
  * [2] K. L. Lueth, _”State of the IoT 2018: Number of IoT devices now at 7B – Market accelerating,”_ in IoT Analytics: Market insight For IoT, 2018 [Online]. Available: https://iot-analytics.com/state-of-the-iot-update-q1-q2-2018-number-of-iot-devices-now-7b/.
  * [3] U. Raza, P. Kulkarni and M. Sooriyabandara, _”Low Power Wide Area Networks: An Overview,”_ in IEEE Commun. Sur. & Tut., vol. 19, 2017.
  * [4] O. Georgiou and U. Raza, _”Low Power Wide Area Network Analysis: Can LoRa Scale?,”_ in IEEE Wireless Commun. Lett., April 2017.
  * [5] ”Low power long range transmitter”, Patent, Jan., 2014. [Online]. Available: https://patents.google.com/patent/US20140219329A1/.
  * [6] LoRa Alliance, _”A technical overview of LoRa and LoRaWAN”,_ [Online]. Available: https://lora-alliance.org/resource-hub/what-lorawantm.
  * [7] T. Elshabrawy and J. Robert, _”Analysis of BER and Coverage Performance of LoRa Modulation under Same Spreading Factor Interference,”_ 2018 IEEE 29th Annu. Int. Sympo. on PIMRC, Bologna, 2018, pp. 1-6.
  * [8] N. E. Rachkidy, A. Guitton and M. Kaneko, _”Decoding Superposed LoRa Signals,”_ 2018 IEEE 43rd Conference on Local Computer Networks (LCN), Chicago, IL, USA, 2018, pp. 184-190.
  * [9] B. Laporte-Fauret, M. A. Ben Temim, G. Ferré, D. Dallet, B. Minger, and L. Fuché, _”An Enhanced LoRa-Like Receiver for the Simultaneous Reception of Two Interfering Signals”_ , 2019 IEEE 30th Annual International Symposium on PIMRC, Istanbul, Turkey, 2019, pp. 1-6.
  * [10] M. A. B. Temim, G. Ferré, B. Laporte-Fauret, D. Dallet, B. Minger and L. Fuché, _”An Enhanced Receiver to Decode Superposed LoRa-like Signals”_ , in IEEE Internet of Things Journal, 2020.
  * [11] U. Noreen, L. Clavier and A. Bounceur, _”LoRa-like CSS-based PHY layer, Capture Effect and Serial Interference Cancellation,”_ European Wireless 2018; 24th European Wireless Conf., Catania, 2018, pp. 1-6.
  * [12] J. G. Andrews, _”Interference cancellation for cellular systems: a contemporary overview,”_ in IEEE Wireless Commun., April 2005.
  * [13] Semtech corporation, ”LoRa Modulation Basics”,, May, 2015. [Online]. Available: https://www.semtech.com.
  * [14] D. Croce, M. Gucciardo, S. Mangione, G. Santaromita and I. Tinnirello, _”Impact of LoRa Imperfect Orthogonality: Analysis of Link-Level Performance,”_ in IEEE Commun. Lett., vol. 22, no. 4, April 2018.