Coverage probability of energy harvesting enabled LoRa networks with stochastic geometry

ABSTRACT

The average transmit power and coverage probability (Pcov) of uplink energy harvesting-enabled long-range networks are investigated in the present paper. Particularly, we model the end-device (EDs) according to the homogeneous Poisson point process while the power beacon is randomly distributed on the circle in the middle of the network. All EDs rely on the harvested energy to perform their operations and transmissions. Under this context, the upper bound of the average transmit power of an end-device is derived in the closed-form expression. The signal-to-noise-ratio condition of the coverage probability is given in the closed-form expression as well. Simulation results are provided to corroborate the accuracy of the derived mathematical framework as well as to feature the impact of some key parameters on the considered metrics. Our findings unveil that increasing either the number of power beacons or their transmit power will monotonically ameliorate the Pcov. Nevertheless, rising the average number of EDs will significantly decline the Pcov's performance.

KEYWORDS:

• Coverage probability
• energy harvesting
• LoRa networks
• performance analysis

1. Introduction

Low power wide area network (LPWAN) has emerged as one of the most effective networks to implement Internet of Things (IoTs) (Raza et al., 2017). More precisely, compared with conventional wireless networks such as cellular networks and local area networks (LAN), that is, WiFi, that require high power consumption as well as well-organized signalling and synchronization. The LPWAN, on the other hand, requires low power consumption at end devices (EDs) and aims for a long transmission distance. It, as a consequence, is suitable for remote terminals like sensors that are hard to connect to the power grid and/or replace the battery. Among all available LPWAN networks, for example, SigFox, LoRa, NB-IoTs, WiFi Halow, etc. (Goursaud & Gorce, 2015). Long-range (LoRa) has gained much attention from both academia and industry. The main advantage of LoRa networks is its patented chirp spread spectrum (CSS) modulation that attains better performance compared with conventional modulations like quadrature amplitude modulation (QAM), and phase shift keying (PSK) (Huan et al., 2022). Another superiority of the LoRa networks is that it allows transmitting over a long transmission distance with ultra-low power consumption by properly configuring a set of parameters such as spreading factor (SF), coding rate (CR), payload, and so on.

Energy harvesting (EH), on the other hand, is one of the most advanced techniques to significantly ameliorate the energy efficiency (EE) of wireless networks (Tu & Di Renzo, 2020). Particularly, by harvesting the ambient radio frequency (RF) signals to recharge the battery, it can not only prolong the devices' lifetime but also be independent of the power grid. Generally, there are two approaches to implementing EH-enabled networks which are power beacon (PB) based and simultaneous wireless information and power transfer (SWIPT) approaches (Lu et al., 2015). The latter permits the receiver to concurrently decode information and absorb energy from the incoming RF signals by properly splitting the received signals in time and/or frequency and/or space. The main drawback of this approach is that instead of directly sending received signals to the decoding circuit like non-EH-enabled networks, they need to part with the harvesting circuit. As a consequence, compared with the performance of the conventional networks in terms of spectral efficiency (SE), the EH-enabled network is not comparable. The former approach, on the other hand, separates the RF signals for harvesting energy and decoding information. More precisely, there are dedicated RF signals broadcast by power beacons to power up wireless devices and are different from conventional RF signals for decoding information. Thus, it does not degrade the performance of the networks like the SWIPT-based approach.

As a result, in the present work, we consider the EH-enabled LoRa networks to facilitate the benefits of both networks and technologies. The performance of LoRa networks was largely investigated in (Bor & Roedig, 2017; Croce et al., 2018; Georgiou & Raza, 2017; Mahmood et al., 2019; Nguyen et al., 2020; Tu et al., 2021, 2022c). More precisely, authors in (Georgiou & Raza, 2017) investigated the coverage probability (Pcov) of an uplink LoRa networks. They, however, ignored the correlation between the signal-to-noise ratios (SNR) and signal-to-interference ratio (SIRs) at the receiver. To overcome such an issue, a novel definition of the Pcov was proposed in (Tu et al., 2022c) where the authors addressed the performance of not only the Pcov but also area spectral efficiency (ASE). A set of LoRa parameters was optimally selected to minimize transmission energy cost in (Bor & Roedig, 2017) while the energy efficiency of the LoRa networks was addressed in (Tu et al., 2021). The authors in (Croce et al., 2018) investigated the imperfect orthogonality between spreading factors. The performance of dual-hop LoRa networks and their scalability were investigated in (Nguyen et al., 2020) and (Mahmood et al., 2019). On the other hand, the performance of the EH-enabled wireless networks was given in (Lakshmi & Jibukumar, 2021; Nasir et al., 2013; Nguyen et al., 2016, 2022, 2023a; Tu & Di Renzo, 2017; Tu et al., 2016). Particularly, Nasir and other authors proposed two protocols to realize the SWIPT-based networks, that is, the time switching and power splitting protocols. The Pcov under the SWIPT mmWave communications was conducted in (Tu & Di Renzo, 2017) while the user selection under the EH-based random networks was investigated in (Nguyen et al., 2016). Nguyen and other authors investigated the trade-off between reliability and security in SWIPT-enabled full-duplex relaying networks (Nguyen et al., 2023a) and a hybrid protocol for SWIPT cooperative networks was proposed in (Lakshmi & Jibukumar, 2021). The receiver diversity under SWIPT cellular networks was investigated in (Tu et al., 2016) while the Pcov of the non-orthogonal multiple access (NOMA) IoT networks with Nakagami-m distribution was derived in (Nguyen et al., 2022).

All the above-mentioned works separately study the performance of the LoRa and EH-enabled networks. The performance of the EH-enabled LoRa networks, on the other hand, is still in the infancy stage with a few works (Hamdi et al., 2022; Hamdi & Qaraqe, 2021). More particular, by optimally allocating the SF and device scheduling, the maximization of the EE in LoRa networks where EDs were powered by a PB was addressed in (Hamdi & Qaraqe, 2021). The power grid minimization was formulated and solved in (Hamdi et al., 2022) by a heuristic algorithm and deep reinforcement learning approach. Regarding the applications of stochastic geometry (SG) and energy harvesting techniques into LoRa networks, only a work from Georgiou and his colleague in (Georgiou et al., 2021) studied such networks. More precisely, they addressed the steady-state distribution of the capacitor voltage and the impact of co-SF interference on the outage probability (OP). We, on the other hand, derived the coverage probability and the average harvested energy. Hence, different from the all above-mentioned works, in the present paper, we investigate the power beacon-based LoRa networks. All power beacon stations are randomly distributed on a circle in the middle of the networks and EDs are modelled according to the homogeneous Poisson point process (PPP). More precisely, the main contributions and novelties are summarized as follows:

• The positions of EDs are modelled according to the Poisson point process while the positions of PBs are randomly distributed on a circle in the middle of the networks.

• We derive the upper bound of the average transmit power of an ED in the closed-form expression. The SNR condition of the Pcov is presented in the closed-form equation as well.

• We employ extensively computer-based simulations to verify the accuracy of the derived mathematical frameworks and to unveil the impact of several key parameters on the considered metrics.

The remainder of the manuscript is organized as follows: Section 2 provides the system model. The main derivations are given in Section 3. Numerical results are provided in Section 4. Section 5 concludes the paper.

2. System model

Let us consider an uplink EH-enabled LoRa network comprises of a single gateway that is located at the centre of the two-dimensional plane with a radius R¹ and multiple end-devices which are distributed around the gateway as shown in Figure 1. The end-devices are modelling according to the homogeneous Poisson point process denoted by ${\mathrm{\Theta }}_{\mathrm{U}}$ with intensity ${\delta }_{\mathrm{U}}$ (Di Renzo et al., 2019b). Additionally, $N\in \mathbb{N}$ power beacons are randomly deployed on the circle with centre at the gateway and radius $\frac{R}{2}$. It is emphasized that the PPP is employed instead of other point processes such as Matérn cluster point processes owing to its mathematical tractability. Moreover, all nodes including the LoRa gateway, power beacon, and end-devices are equipped with only a single antenna. The considered network is operated at the industrial, scientific, and medical (ISM) band. Yet the interference from other technologies operating with the same band is not taken into consideration owing to the advantages of CSS modulations. Different from the conventional LoRaWAN protocol, we adopt the slotted-ALOHA protocol at the medium access control (MAC) layer since its performance is better than the conventional one (Beltramelli et al., 2021).

Figure 1. The considered power beacon-based LoRa networks.

2.1. Channel modelling

In this work, the transmission link from EDs to the gateway is suffered from both small-scale fading and large-scale path loss. The impact of the shadowing can be implicitly studied by modifying the density of the EDs (Di Renzo et al., 2018).

2.1.1. Small-scale fading

Considering a generic transmission link between EDs and the gateway, the channel coefficient denoted by h is modelled by a complex Gaussian random variable (RV) with zero mean and Ω variance, that is, $h\in \mathcal{CN}(0,\mathrm{\Omega })$. As a consequence, the channel gain denoted by $|h{|}^2$ is followed by an exponential distribution with mean Ω. We further assume that the small-scale fading is a block flat fading thus, it remains stable for a whole time slot and changes independently between time slots (Tu et al., 2023b).

2.1.2. Large-scale path-loss

Besides the effects of the small-scale fading, the transmission link from EDs to the gateway is also subjected to the large-scale path-loss denoted by Λ and is calculated as follows: (1) $\mathrm{\Lambda }={\mathcal{K}}_{0}max\{1,d^{\alpha }\},$(1) where $max\{.\}$ is the maximum function and ${\mathcal{K}}_0=(\frac{4\pi }{\iota }{)}^2$ is the path-loss constant, $\alpha \in (2,8)$ is the path-loss exponent relying on the transmission environment. For example, under the urban environment α will be generally between 3.5−4.5 while in the rural areas, its value will approach 2. $\iota =\frac{c}{f_c}$ is the wavelength, $f_c$ (in Hz) and $c=3\times {10}^8$ are the carrier frequency and the speed of light. d is the transmission distance from EDs to the gateway.

Remark 2.1

It is certain that the adopted large-scale path-loss model has overcome the singularity issue in the popular unbounded path-loss model where the transmission distance approaches zero, and the received power then approaches infinity (Aravanis et al., 2019).

2.2. Spreading factor allocation

In this work, the fair-collision probability allocation is adopted. Under this scheme, the probability that an ED is allocated to the kth SF, $k\in \{7,\dots ,12\}$, is computed by the following equation: (2) $p_k=\frac{k/2^k}{\sum _{l=7}^{12}l/2^l},k\in \{7,\dots ,12\}.$(2)

Remark 2.2

It is emphasized that the adopted fair-collision probability SF allocation achieves better performance compared with the popular distance-based approach (Reynders et al., 2017).

2.3. Energy harvesting modelling

We assume that the ED is not connected to the power grid and counts solely on the harvested energy from all power beacons during sleep mode. Let us denote ${\mathcal{E}}_{k,o}$ by the harvested energy of the oth ED utilized SFk, it is formulated as follows (Nguyen et al., 2017): (3) ${\mathcal{E}}_{k,o}=\sum _{i=1}^N\frac{\chi P_P{\left|h_{i,o}\right|}^2{T}_{\mathrm{sleep},k}}{{\mathrm{\Lambda }}_{i,o}},$(3) where $P_P$ is the transmit power of the power beacon. Here, we assume that all power beacons emit the same transmit power. ${\left|h_{i,o}\right|}^2$ is the channel gain from the ith power beacon to the oth ED and $T_{\mathrm{sleep},k}=T_{\mathrm{all},k}-T_{\mathrm{trans},k}$ is the sleeping period of the ED utilized SFk and is the function of the transmission duration $T_{\mathrm{trans},k}$ and the total transmission duration $T_{\mathrm{all},k}=max\left\{T_{\mathrm{not}},\frac{1}{\rho }{T}_{\mathrm{trans},k}\right\}$. Here, $T_{\mathrm{not}}$ is the notification duration between two transmissions, and ρ is the duty cycle. Generally, ρ is defined according to government regulation and is not greater than 1%. $\chi \in \left(0,1\right)$ is the energy conversion efficiency.

2.4. Signal-to-Noise ratio and signal-to-interference ratio

2.4.1. Signal-to-Noise ratio

The signal-to-noise ratio at the gateway from the desired ED denoted by subscript ‘0’ operated at SFk is formulated as (4) $\begin{array}{rl}{\mathrm{SNR}}_k& =\frac{P_{0,k}{\left|h_0\right|}^2}{{\mathrm{\Lambda }}_0{\sigma }^2},\\ P_{0,k}& =\frac{{\mathcal{E}}_{k,o}}{T_{\mathrm{trans},k}}.\end{array}$(4) Here $P_{0,k}$ is the transmit power of the intended ED and is a function of the harvested energy. It is noted that we skip the energy consumption for other activities since it is negligible compared with the transmit power. ${\sigma }^2=-174+\mathrm{NF}+10{\log }_{10}\left\{\mathrm{Bw}\right\}$ is the noise variance of the additive white Gaussian noise (AWGN); NF (in dB) is the noise figure of the gateway, BW (in Hz) is the transmission bandwidth, $\log \{.\}$ is the logarithm function; ${\left|h_0\right|}^2$ is the channel gain from desired ED to the gateway.

2.4.2. Signal-to-Interference Ratio

Besides the SNR, the received signals at the gateway also suffer from interference from co-SF interference. The signal-to-interference ratio is formulated as follows: (5) ${\mathrm{SIR}}_k=\frac{P_{0,k}{\left|h_0\right|}^2/{\mathrm{\Lambda }}_0}{\sum _{i\in {\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A}{P}_{i,k}{\left|h_i\right|}^2/{\mathrm{\Lambda }}_{i,k}},$(5) where ${\left|h_i\right|}^2$ is the channel gain from the ith ED utilized SFk, $P_{i,k}=\frac{{\mathcal{E}}_{k,i}}{T_{\mathrm{trans},k}}$ is its transmit power, and ${\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A=\rho p_k{\delta }_{\mathrm{U}}$ is the set of active EDs utilized SFk.

2.5. Coverage probability

In this work, we study the coverage probability², Pcov, of the ED. The coverage probability is the probability that the received signals at the gateway are greater than a predefined threshold. It is noted that different from other wireless networks the predefined threshold of LoRa networks is a discrete set and relies on the utilized SF. Additionally, the Pcov comprises two conditions, the SNR and SIR conditions, instead of a single condition in the form of the signal-to-interference-plus-noise (SINR). Mathematical speaking, the Pcov under SFk is formulated as follows (Di Renzo et al., 2019a): (6) $P_{\mathrm{cov}}\left(k\right)=Pr\left\{{\mathrm{SNR}}_{\mathrm{k}}\ge q_k,{\mathrm{SIR}}_{\mathrm{k}}\ge {ϵ}_k\right\}.$(6) Here $Pr\{.\}$ is the probability operator; $q\in \left\{-6,-9,-12,-15,-17.5,-20\right\}$ (in dBm), $k\in \{7,\dots ,12\}$ is the target threshold; ${ϵ}_k=ϵ=1\mathrm{dB},\mathrm{\forall k}$ is the co-SF interference threshold and is equal to 1 dB (Croce et al., 2018, Table II).

3. Performance analysis

Before going to compute the Pcov, the two following useful Lemmas that address the performance of the average transmit power of each ED and the characteristics of the ratio of the channel gain and path-loss from an ED to the gateway are given as follows:

Lemma 3.1

Let us denote ${\mathcal{P}}_k$ by the upper bound of the average transmit power from (3(3) ${\mathcal{E}}_{k,o}=\sum _{i=1}^N\frac{\chi P_P{\left|h_{i,o}\right|}^2{T}_{\mathrm{sleep},k}}{{\mathrm{\Lambda }}_{i,o}},$(3) ) and (4(4) $\begin{array}{rl}{\mathrm{SNR}}_k& =\frac{P_{0,k}{\left|h_0\right|}^2}{{\mathrm{\Lambda }}_0{\sigma }^2},\\ P_{0,k}& =\frac{{\mathcal{E}}_{k,o}}{T_{\mathrm{trans},k}}.\end{array}$(4) ) of ED utilized SFk, it then computes as (7) ${\mathcal{P}}_k=N\frac{\chi P_P{T}_{\mathrm{sleep},k}}{T_{\mathrm{trans},k}}\frac{2\mathrm{\Omega }}{R}\left(1+\frac{1}{\alpha -1}\left(1-{\left(\frac{R}{2}\right)}^{1-a}\right)\right),$(7)

where $\mathbb{E}\{.\}$ is the expectation operator.

Proof.

The proof is given in Appendix A.

Lemma 3.2

Let us denote random variable $Y=\frac{{\left|h\right|}^2}{max\left\{1,d^{\alpha }\right\}}$, its complementary cumulative distribution function (CCDF), CDF, and PDF are computed as follows: (8) $\begin{array}{rl}{\overline{F}}_Y\left(y\right)& =\frac{1}{R^2}\exp \left(-\frac{y}{\mathrm{\Omega }}\right)+\frac{1}{\alpha }\frac{2}{R^2}{\left(\frac{y}{\mathrm{\Omega }}\right)}^{-\frac{2}{\alpha }}\left[\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{y}{\mathrm{\Omega }}\right)-\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{y{R}^{\alpha }}{\mathrm{\Omega }}\right)\right]\\ F_Y\left(y\right)& =1-\frac{1}{R^2}\exp \left(-\frac{y}{\mathrm{\Omega }}\right)-\frac{1}{\alpha }\frac{2}{R^2}{\left(\frac{y}{\mathrm{\Omega }}\right)}^{-\frac{2}{\alpha }}\left[\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{y}{\mathrm{\Omega }}\right)-\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{y{R}^{\alpha }}{\mathrm{\Omega }}\right)\right]\\ f_Y\left(y\right)& =\frac{1}{\mathrm{\Omega }{R}^2}\exp \left(-\frac{y}{\mathrm{\Omega }}\right)+\frac{4}{{\alpha }^2{R}^2}{\left(\frac{1}{\mathrm{\Omega }}\right)}^{-\frac{2}{\alpha }}\left(y^{-1-\frac{2}{\alpha }}\right)\\ & \times \left[\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{y}{\mathrm{\Omega }}\right)-\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{y{R}^{\alpha }}{\mathrm{\Omega }}\right)\right]+\frac{2}{\alpha R^2}\frac{1}{y}\\ & \times \left[\exp \left(-\frac{y}{\mathrm{\Omega }}\right)-R^2\exp \left(-\frac{y{R}^{\alpha }}{\mathrm{\Omega }}\right)\right].\end{array}$(8) Here ${\overline{F}}_Y\left(y\right)$ is the CCDF of RV Y.

Proof.

The proof is given in Appendix B.

Remark 3.1

It is emphasized that the outcomes in Lemma 3.2 can be straightforwardly applied to the Nakagami-m distribution provided that m is an integer number, that is, $m\in \mathbb{N}$.

Having obtained the average transmit power and the CCDF, CDF, and PDF of the ratio of the channel gain and large-scale path-loss, we are ready to compute the considered metric, that is, the coverage probability. Let us begin with the definition of Pcov in (6(6) $P_{\mathrm{cov}}\left(k\right)=Pr\left\{{\mathrm{SNR}}_{\mathrm{k}}\ge q_k,{\mathrm{SIR}}_{\mathrm{k}}\ge {ϵ}_k\right\}.$(6) ) as follows: (9) $\begin{array}{rl}{P}_{\mathrm{cov}}\left(k\right)& =Pr\left\{{\mathrm{SNR}}_{\mathrm{k}}\ge q_k,{\mathrm{SIR}}_{\mathrm{k}}\ge {ϵ}_k\right\}\\ & =Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2}{{\mathrm{\Lambda }}_0{\sigma }^2}\ge q_k,\frac{P_{0,k}{\left|h_0\right|}^2/{\mathrm{\Lambda }}_0}{\sum _{i\in {\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A}{P}_{i,k}{\left|h_i\right|}^2/{\mathrm{\Lambda }}_{i,k}}\ge {ϵ}_k\right\}\\ & \stackrel{\left(a\right)}{=}Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2}{{\mathrm{\Lambda }}_0{\sigma }^2}\ge q_k\right\}Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2/{\mathrm{\Lambda }}_0}{\sum _{i\in {\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A}{P}_{i,k}{\left|h_i\right|}^2/{\mathrm{\Lambda }}_{i,k}}\ge {ϵ}_k\right\}\\ & ={\mathcal{I}}_1{\mathcal{I}}_2,\end{array}$(9) where $\left(a\right)$ is held by assuming that the SNR and SIR are uncorrelated random variables³. The integration ${\mathcal{I}}_1$ and ${\mathcal{I}}_2$ refer to the SNR and SIR conditions of the Pcov and are computed as follows: (10) $\begin{array}{rl}{\mathcal{I}}_1& =Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2}{{\mathrm{\Lambda }}_0{\sigma }^2}\ge q_k\right\}\approx Pr\left\{\frac{{\mathcal{P}}_k}{{\mathcal{K}}_0{\sigma }^2}{Y}_0\ge q_k\right\}={\overline{F}}_Y\left(\frac{{\mathcal{K}}_0{q}_k{\sigma }^2}{{\mathcal{P}}_k}\right)\\ & \stackrel{\left(a\right)}{=}\exp \left(-\frac{{\mathcal{K}}_0{q}_k{\sigma }^2}{{\mathcal{P}}_k}\right)\frac{1}{R^2}+\frac{1}{\alpha }\frac{2}{R^2}{\left(\frac{{\mathcal{K}}_0{q}_k{\sigma }^2}{{\mathcal{P}}_k}\right)}^{-\frac{2}{\alpha }}\\ & \times \left[\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{{\mathcal{K}}_0{q}_k{\sigma }^2}{{\mathcal{P}}_k}\right)-\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{{\mathcal{K}}_0{q}_k{\sigma }^2}{{\mathcal{P}}_k}{R}^{\alpha }\right)\right],\end{array}$(10) where $Y_0=\frac{{\left|h_0\right|}^2}{max\left\{1,d_0^{\alpha }\right\}}$ and $\left(a\right)$ is achieved with the help of Lemma 3.2 and ${\mathcal{P}}_k$ comes from Lemma 3.1. The second integration, ${\mathcal{I}}_2$ is computed as follows: (11) $\begin{array}{rl}{\mathcal{I}}_2& =Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2/{\mathrm{\Lambda }}_0}{\sum _{i\in {\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A}{P}_{i,k}{\left|h_i\right|}^2/{\mathrm{\Lambda }}_{i,k}}\ge {ϵ}_k\right\}\stackrel{\left(a\right)}{\le }Pr\left\{\frac{{\mathcal{P}}_k{\left|h_0\right|}^2/{\mathrm{\Lambda }}_0}{{\mathcal{P}}_k\underset{i\in {\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A}{max}\left\{{\left|h_i\right|}^2/{\mathrm{\Lambda }}_{i,k}\right\}}\ge {ϵ}_k\right\}\\ & =Pr\left\{\frac{{\mathcal{P}}_k{Y}_0}{{\mathcal{P}}_k\underset{i\in {\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A}{max}\left\{Y_i\right\}}\ge {ϵ}_k\right\}=Pr\{\frac{Y_0}{Z}\ge {ϵ}_k\}\\ & ={\int }_0^{\mathrm{\infty }}{\overline{F}}_{Y_0}\left(z{ϵ}_k\right)f_Z\left(z\right)\mathrm{d}z={\int }_0^{\mathrm{\infty }}[\frac{1}{R^2}\exp (-\frac{z{ϵ}_k}{\mathrm{\Omega }})+\frac{1}{\alpha }\frac{2}{R^2}\\ & \times {(-\frac{z{ϵ}_k}{\mathrm{\Omega }})}^{-\frac{2}{\alpha }}\left[\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{z{ϵ}_k}{\mathrm{\Omega }}\right)-\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{z{ϵ}_k{R}^{\alpha }}{\mathrm{\Omega }}\right)\right]]\\ & (\times \sum _{i=1}^O(\frac{1}{{\mathrm{\Omega }}_i{R}^2}\exp \left(-\frac{z}{{\mathrm{\Omega }}_i}\right)+\frac{4}{{\alpha }^2{R}^2}{\left(\frac{1}{{\mathrm{\Omega }}_i}\right)}^{-\frac{2}{\alpha }}\left(z^{-1-\frac{2}{\alpha }}\right)[\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{y}{{\mathrm{\Omega }}_i}\right)\\ & -\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{y{R}^{\alpha }}{{\mathrm{\Omega }}_i}\right)]+\frac{2}{\alpha R^2}\left[\frac{1}{z}\exp \left(-\frac{z}{{\mathrm{\Omega }}_i}\right)-{\left(R^{\alpha }\right)}^{\frac{2}{\alpha }}\exp \left(-\frac{z{R}^{\alpha }}{{\mathrm{\Omega }}_i}\right)\right])\mathrm{d}z\\ & \prod _{n=1,n\ne i}^O\left(1-\frac{1}{R^2}\exp \left(-\frac{z}{{\mathrm{\Omega }}_n}\right)-\frac{1}{\alpha }\frac{2}{R^2}{\left(\frac{z}{{\mathrm{\Omega }}_n}\right)}^{-\frac{2}{\alpha }}\left[\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{z}{{\mathrm{\Omega }}_n}\right)-\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{z{R}^{\alpha }}{{\mathrm{\Omega }}_n}\right)\right]\right)],\end{array}$(11) where $\left(a\right)$ is held by approximating the aggregate interference by the strongest interferer, $Z=\underset{i\in {\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A}{max}\left\{Y_i\right\}$, O is the number of active EDs. The CDF and PDF of RV Z are calculated by $F_Z\left(z\right)=\prod _{i=1}^O{F}_{Y_i}\left(y\right)$ and $f_Z\left(z\right)=\sum _{i=1}^O{f}_{Y_i}\left(y\right)\prod _{n=1,n\ne i}^O{F}_{Y_n}\left(y\right)$. The last equation can be effectively computed by employing popular software such as Matlab, Mathematica, and Python.

4. Numerical results

In this section, computer-based simulation results are provided to justify the precision of the developed mathematical framework as well as to reveal the influences of some key parameters such as the average number of end devices, the number of power beacons, and the transmit power of the PB. Without loss of generality, the following set of parameters are used in this section: Bw = 125 kHz, α = 2.5, $f_c$ = 920 MHz, $P_P$ = 30 dBm, N = 5, NF = 6 dB, R = 2000 m, Ω = 3.5, ${\delta }_U=\frac{\overline{N}}{\pi R^2}$, $\overline{N}$ = 500, ρ = 0.01, $\chi =0.55$, $T_{\mathrm{trans},k}=3\mathrm{seconds}\mathrm{\forall }k$, $T_{\mathrm{not},k}=1800\mathrm{seconds}\mathrm{\forall }k$, respectively.

Figure 2 stretches the average transmit power of the ED w.r.t. the number of power beacons. We first observe that there is a good agreement between the derived mathematical framework and the simulation results. Second, the average transmit power is a monotonic increasing function of the N which can be straightforwardly seen in (7(7) ${\mathcal{P}}_k=N\frac{\chi P_P{T}_{\mathrm{sleep},k}}{T_{\mathrm{trans},k}}\frac{2\mathrm{\Omega }}{R}\left(1+\frac{1}{\alpha -1}\left(1-{\left(\frac{R}{2}\right)}^{1-a}\right)\right),$(7) ). Third, the higher the path-loss exponent the smaller the average transmit power. This phenomenon can be explained that with the larger the path-loss exponent the higher the path-loss and the smaller the receiver power at the ED.

Figure 2. Average transmit power vs. number of power beacons, N with various values of path-loss exponent. Solid lines are computed by (7(7) ${\mathcal{P}}_k=N\frac{\chi P_P{T}_{\mathrm{sleep},k}}{T_{\mathrm{trans},k}}\frac{2\mathrm{\Omega }}{R}\left(1+\frac{1}{\alpha -1}\left(1-{\left(\frac{R}{2}\right)}^{1-a}\right)\right),$(7) ) while markers are from Monte-Carlo simulations.

Figure 3 illustrates the CDF of RV X defined in (A1(A1) $\begin{array}{rl}\mathbb{E}\left\{P_{i,k}\right\}& =\sum _{i=1}^N\frac{\chi P_P{\left|h_{i,o}\right|}^2{T}_{\mathrm{sleep},k}}{max\left\{1,{\Vert {\mathbf{p}}_o-{\mathbf{p}}_i\Vert }^{\alpha }\right\}{T}_{\mathrm{trans},k}}\le {\mathcal{P}}_k=\sum _{i=1}^N\frac{\chi P_P{\left|h_{i,o}\right|}^2{T}_{\mathrm{sleep},k}}{max\left\{1,t_o^{\alpha }\right\}{T}_{\mathrm{trans},k}}\\ & =N\frac{\chi P_P{T}_{\mathrm{sleep},k}}{T_{\mathrm{trans},k}}\mathbb{E}\left\{X=\frac{{\left|h_{i,o}\right|}^2}{max\left\{1,t_o^{\alpha }\right\}}\right\}.\end{array}$(A1) ) with different values of the network radius R. We see that the larger the network radius the smaller the CDF. The rationale behind this observation is that the transmission distance is getting further hence the RV X becomes smaller so achieves a higher probability. We experience again that the curves from the mathematical framework absolutely align with the Monte Carlo results.

Figure 3. CDF of RV X in (A1(A1) $\begin{array}{rl}\mathbb{E}\left\{P_{i,k}\right\}& =\sum _{i=1}^N\frac{\chi P_P{\left|h_{i,o}\right|}^2{T}_{\mathrm{sleep},k}}{max\left\{1,{\Vert {\mathbf{p}}_o-{\mathbf{p}}_i\Vert }^{\alpha }\right\}{T}_{\mathrm{trans},k}}\le {\mathcal{P}}_k=\sum _{i=1}^N\frac{\chi P_P{\left|h_{i,o}\right|}^2{T}_{\mathrm{sleep},k}}{max\left\{1,t_o^{\alpha }\right\}{T}_{\mathrm{trans},k}}\\ & =N\frac{\chi P_P{T}_{\mathrm{sleep},k}}{T_{\mathrm{trans},k}}\mathbb{E}\left\{X=\frac{{\left|h_{i,o}\right|}^2}{max\left\{1,t_o^{\alpha }\right\}}\right\}.\end{array}$(A1) ) with various values of the network radius R. Solid lines are computed by (A2(A2) $\begin{array}{rl}{F}_X\left(x\right)& =Pr\left\{X=\frac{{\left|h_{i,o}\right|}^2}{max\left\{1,t_o^{\alpha }\right\}}\le x\right\}\stackrel{\left(a\right)}{=}{\int }_0^{R/2}\left(1-\exp \left(-\frac{xmax\left\{1,t^{\alpha }\right\}}{\mathrm{\Omega }}\right)\right)\frac{2}{R}\mathrm{d}t\\ & =1-\frac{2}{R}{\int }_0^{R/2}\exp \left(-\frac{xmax\left\{1,t^{\alpha }\right\}}{\mathrm{\Omega }}\right)\mathrm{d}t\\ & =1-\frac{2}{R}\left[{\int }_0^1\exp \left(-\frac{x}{\mathrm{\Omega }}\right)\mathrm{d}t+{\int }_1^{R/2}\exp \left(-\frac{x{t}^{\alpha }}{\mathrm{\Omega }}\right)\mathrm{d}t\right]\\ & \stackrel{\left(b\right)}{=}1-\frac{2}{R}\exp \left(-\frac{x}{\mathrm{\Omega }}\right)-\frac{2}{R}\frac{1}{\alpha }{\left(\frac{\mathrm{\Omega }}{x}\right)}^{\frac{1}{\alpha }}\left[\mathrm{\Gamma }\left(\frac{1}{\alpha },\frac{x}{\mathrm{\Omega }}\right)-\mathrm{\Gamma }\left(\frac{1}{\alpha },\frac{x}{\mathrm{\Omega }}{\left(\frac{R}{2}\right)}^{\alpha }\right)\right],\end{array}$(A2) ) while markers are from Monte-Carlo simulations.

The influences of the transmit power of ED on the performance of coverage probability is given in Figure 4. We discover that the mathematical framework denoted by ‘Framework’ is a lower bound of both simulation results. The curves with dot red markers are simulated results that consider the strongest interference while the curves with diamond black markers consider the aggregate interference. All curves plotted by Monte Carlo simulations take into account the correlation between the SNR and SIR in (9(9) $\begin{array}{rl}{P}_{\mathrm{cov}}\left(k\right)& =Pr\left\{{\mathrm{SNR}}_{\mathrm{k}}\ge q_k,{\mathrm{SIR}}_{\mathrm{k}}\ge {ϵ}_k\right\}\\ & =Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2}{{\mathrm{\Lambda }}_0{\sigma }^2}\ge q_k,\frac{P_{0,k}{\left|h_0\right|}^2/{\mathrm{\Lambda }}_0}{\sum _{i\in {\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A}{P}_{i,k}{\left|h_i\right|}^2/{\mathrm{\Lambda }}_{i,k}}\ge {ϵ}_k\right\}\\ & \stackrel{\left(a\right)}{=}Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2}{{\mathrm{\Lambda }}_0{\sigma }^2}\ge q_k\right\}Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2/{\mathrm{\Lambda }}_0}{\sum _{i\in {\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A}{P}_{i,k}{\left|h_i\right|}^2/{\mathrm{\Lambda }}_{i,k}}\ge {ϵ}_k\right\}\\ & ={\mathcal{I}}_1{\mathcal{I}}_2,\end{array}$(9) ). It is no doubt that the curves considering the strongest interferer are an upper bound of the aggregate interference. Nevertheless, the gap between these cases is minor. Thus, it is reasonable to approximate the aggregate interference by the strongest interferer in LoRa networks. From Figure 4, the Pcov of SF12 achieves a higher performance compared with the ones utilized by SF7 under the current setup. It is expected that increasing the transmit power of the power beacon is beneficial for the Pcov. However, there exists an upper limit of the Pcov if $P_P$ keeps increasing.

Figure 4. Coverage probability vs. transmit power of power beacon, $P_P$. Solid lines are computed by (9(9) $\begin{array}{rl}{P}_{\mathrm{cov}}\left(k\right)& =Pr\left\{{\mathrm{SNR}}_{\mathrm{k}}\ge q_k,{\mathrm{SIR}}_{\mathrm{k}}\ge {ϵ}_k\right\}\\ & =Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2}{{\mathrm{\Lambda }}_0{\sigma }^2}\ge q_k,\frac{P_{0,k}{\left|h_0\right|}^2/{\mathrm{\Lambda }}_0}{\sum _{i\in {\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A}{P}_{i,k}{\left|h_i\right|}^2/{\mathrm{\Lambda }}_{i,k}}\ge {ϵ}_k\right\}\\ & \stackrel{\left(a\right)}{=}Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2}{{\mathrm{\Lambda }}_0{\sigma }^2}\ge q_k\right\}Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2/{\mathrm{\Lambda }}_0}{\sum _{i\in {\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A}{P}_{i,k}{\left|h_i\right|}^2/{\mathrm{\Lambda }}_{i,k}}\ge {ϵ}_k\right\}\\ & ={\mathcal{I}}_1{\mathcal{I}}_2,\end{array}$(9) ) while markers are from Monte-Carlo simulations. Dot red markers consider the strongest interferer while diamond black markers consider the aggregate interference. All simulation curves take into account the correlation between SNR and SIR.

Figure 5 addresses the performance of the Pcov w.r.t the number of power beacons. Inspecting this figure, we observe that the Pcov behaviours w.r.t. the number of power beacon N is similar to $P_P$ that it is a monotonic increasing function. Nonetheless, the increasing pace is nonlinear. More precisely, the Pcov dramatically goes up when N is small, that is, N<10, it then rises at a lower pace. We also witness a negligible gap between curves considering all interference, the strongest interferer, and the framework, especially for SF12.

Figure 5. Coverage probability vs. the number of power beacon. Solid lines are computed by (9(9) $\begin{array}{rl}{P}_{\mathrm{cov}}\left(k\right)& =Pr\left\{{\mathrm{SNR}}_{\mathrm{k}}\ge q_k,{\mathrm{SIR}}_{\mathrm{k}}\ge {ϵ}_k\right\}\\ & =Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2}{{\mathrm{\Lambda }}_0{\sigma }^2}\ge q_k,\frac{P_{0,k}{\left|h_0\right|}^2/{\mathrm{\Lambda }}_0}{\sum _{i\in {\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A}{P}_{i,k}{\left|h_i\right|}^2/{\mathrm{\Lambda }}_{i,k}}\ge {ϵ}_k\right\}\\ & \stackrel{\left(a\right)}{=}Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2}{{\mathrm{\Lambda }}_0{\sigma }^2}\ge q_k\right\}Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2/{\mathrm{\Lambda }}_0}{\sum _{i\in {\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A}{P}_{i,k}{\left|h_i\right|}^2/{\mathrm{\Lambda }}_{i,k}}\ge {ϵ}_k\right\}\\ & ={\mathcal{I}}_1{\mathcal{I}}_2,\end{array}$(9) ) while markers are from Monte-Carlo simulations. Dot red markers consider the strongest interferer while diamond black markers consider the aggregate interference. All simulation curves take into account the correlation between SNR and SIR.

The impact of the inter-arrival between two transmissions on the Pcov performance is provided in Figure 6. It is evident that the longer the notification between two transmissions (or the fewer the number of transmissions), the better the Pcov. It can be effortlessly explained that the ED will harvest more energy when the duration between two transmissions extends, boosting the transmit power and the Pcov. Moreover, it is interesting to point out that the Pcov changes from a noise-limited regime to an interference-limited regime when $T_{\mathrm{not}}$ increases. Particularly, when $T_{\mathrm{not}}$ is small, that is, $T_{\mathrm{not}}<3$ minutes, the Pcov is in a noise-limited regime because the amount of harvested energy is modest which leads to a limited transmit power. Nonetheless, when $T_{\mathrm{not}}$ is getting larger, the system gradually changes to an interference-limited regime. Because the harvested energy is enough to overcome the negative impact of the background noise. Finally, we experience again that the larger the SF the higher the Pcov.

Figure 6. Coverage probability vs. the notification duration between two transmissions. Solid lines are computed by (9(9) $\begin{array}{rl}{P}_{\mathrm{cov}}\left(k\right)& =Pr\left\{{\mathrm{SNR}}_{\mathrm{k}}\ge q_k,{\mathrm{SIR}}_{\mathrm{k}}\ge {ϵ}_k\right\}\\ & =Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2}{{\mathrm{\Lambda }}_0{\sigma }^2}\ge q_k,\frac{P_{0,k}{\left|h_0\right|}^2/{\mathrm{\Lambda }}_0}{\sum _{i\in {\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A}{P}_{i,k}{\left|h_i\right|}^2/{\mathrm{\Lambda }}_{i,k}}\ge {ϵ}_k\right\}\\ & \stackrel{\left(a\right)}{=}Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2}{{\mathrm{\Lambda }}_0{\sigma }^2}\ge q_k\right\}Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2/{\mathrm{\Lambda }}_0}{\sum _{i\in {\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A}{P}_{i,k}{\left|h_i\right|}^2/{\mathrm{\Lambda }}_{i,k}}\ge {ϵ}_k\right\}\\ & ={\mathcal{I}}_1{\mathcal{I}}_2,\end{array}$(9) ) while markers are from Monte-Carlo simulations. Dot red markers consider the strongest interferer while diamond black markers consider the aggregate interference. All simulation curves take into account the correlation between SNR and SIR.

Figure 7 studies the Pcov performance regarding the average number of EDs. We see that increasing $\overline{N}$ will decline the Pcov as the interference scales up with $\overline{N}$. The gap between the proposed mathematical framework and the case considering the maximal interference is indistinguishable for both SF7 and SF11.

Figure 7. Coverage probability vs. the average number of EDs, that is, ${\overline{N}}_U={\delta }_U\left(\pi R^2\right)$. Solid lines are computed by (9(9) $\begin{array}{rl}{P}_{\mathrm{cov}}\left(k\right)& =Pr\left\{{\mathrm{SNR}}_{\mathrm{k}}\ge q_k,{\mathrm{SIR}}_{\mathrm{k}}\ge {ϵ}_k\right\}\\ & =Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2}{{\mathrm{\Lambda }}_0{\sigma }^2}\ge q_k,\frac{P_{0,k}{\left|h_0\right|}^2/{\mathrm{\Lambda }}_0}{\sum _{i\in {\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A}{P}_{i,k}{\left|h_i\right|}^2/{\mathrm{\Lambda }}_{i,k}}\ge {ϵ}_k\right\}\\ & \stackrel{\left(a\right)}{=}Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2}{{\mathrm{\Lambda }}_0{\sigma }^2}\ge q_k\right\}Pr\left\{\frac{P_{0,k}{\left|h_0\right|}^2/{\mathrm{\Lambda }}_0}{\sum _{i\in {\mathrm{\Theta }}_{U,k\mathrm{\setminus }0}^A}{P}_{i,k}{\left|h_i\right|}^2/{\mathrm{\Lambda }}_{i,k}}\ge {ϵ}_k\right\}\\ & ={\mathcal{I}}_1{\mathcal{I}}_2,\end{array}$(9) ) while markers are from Monte-Carlo simulations. Dot red markers consider the strongest interferer while diamond black markers consider the aggregate interference. All simulation curves take into account the correlation between SNR and SIR.

5. Conclusion

We studied the performance of the average transmit power and the coverage probability of an uplink EH-enabled LoRa network. More particular, we presented in explicit form the average transmit power and the SNR condition of the Pcov. Numerical results revealed that increasing the number of power beacons and their transmit power is beneficial for the Pcov. Nonetheless, rising the average number of EDs is harmful for this metric. The present work can be extended in several directions. One of the feasible directions is to apply the Fountain codes to enhance the transmission efficiency (Ngoc et al., 2022; Tu et al., 2022d, 2023a). Another interesting direction is to apply MIMO and/or transmit/receive diversity techniques to significantly ameliorate the system performance (Hossain, 2022; Mahey & Malhotra, 2015). Other possible extensions can be considered such as mmWave communications to scale up the spectral efficiency, reconfigurable intelligent surfaces (RIS) (Al Hajj et al., 2020; Chien et al., 2022; Tu et al., 2015; Van Chien et al., 2021). Applying deep learning (DL) into LoRa networks is also an interesting extension (Tu et al., 2022a).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research is funded by Posts and Telecommunications Institute of Technology [grant number 19-2023-HV-TH2].

Notes on contributors

Thi-Tuyet-Hai Nguyen

Thi-Tuyet-Hai Nguyen (Email: [email protected]) got the Ph.D. degree from La Rochelle University, La Rochelle, France, in 2020. She has been working at the Faculty of Information Technology 2, Posts and Telecommunications Institute of Technology, HCM campus, Vietnam since 2012. Her research interests are natural language processing, machine learning, information retrieval, chatbot, and LoRa networks.

Tan N. Nguyen

Tan N. Nguyen (Email: [email protected]) was born in 1986 in Nha Trang City, Vietnam. He received a BS degree in electronics in 2008 from Ho Chi Minh University of Natural Sciences and an MS degree in telecommunications engineering in 2012 from Vietnam National University. He received a Ph.D. in communications technologies in 2019 from the Faculty of Electrical Engineering and Computer Science at VSB – Technical University of Ostrava, Czech Republic. He joined the Faculty of Electrical and Electronics Engineering of Ton Duc Thang University, Vietnam, in 2013, and since then has been lecturing. His major interests are cooperative communications, cognitive radio, physical layer security and signal processing.

Tran Trung Duy

Tran Trung Duy (Email: [email protected]) received the Ph.D degree in electrical engineering from University of Ulsan, South Korea in 2013. In 2013, he joined Posts and Telecommunications Institute of Technology, HoChiMinh city campus (PTIT-HCM), as a lecturer. From 2022, he is an associate Professor of Wireless Communications at PTIT-HCM. His major research interests are cooperative communications, cooperative multi-hop, cognitive radio, physical-layer security, energy harvesting, hardware impairments and Fountain codes.

Nguyen Hong Son

Nguyen Hong Son (Email: [email protected]) received his B.Sc. in Computer Engineering from The University of Technology in HCM city, and his M.Sc. and Ph.D. in Communication Engineering from the Post and Telecommunication Institute of Technology in Hanoi, Vietnam. His research interests include communication engineering, machine learning, network security, and cloud computing.

Tan Hanh

Tan Hanh (Email: [email protected]) received the PhD degree from Grenoble Institute of Technology, France. Currently, he is vice president of Posts and Telecommunications Institute of Technology (PTIT), in charge of Ho Chi Minh city campus. His research interests are machine learning, information retrieval, and data mining.

Minh Bui Vu

Minh Bui Vu (Email: [email protected]) was born on March 02, 1991 in Dong Nai, Vietnam. He graduated in Electrical and Electronic Engineering in 2014 from Nguyen Tat Thanh University, Ho Chi Minh City, Vietnam. End of 2014, he joined the Faculty of Electrical and Electronics of Nguyen Tat Thanh University as Laboratory-Practice management, until in 2017 he was a lecturer. In 2019, he received Master degree in Engineering Electrical from Ho Chi Minh City University of Technology and Education. His major research interests are Artificial Neuron Network, Fuzzy Logic, Wireless Networks.

Lam-Thanh Tu

Lam-Thanh Tu (Email: [email protected]) received the Ph.D. degree from the University of Paris Sud, Paris-Saclay University, France, in 2018. From 2022, he has been with the Faculty of Electrical and Electronics Engineering, at Ton Duc Thang University, Vietnam. His research interests include stochastic geometry, LoRa networks, reconfigurable intelligent surfaces, covert communications, and artificial intelligence applications for wireless communications.

Notes

1 Here we consider a finite network radius $R<\mathrm{\infty }$, if the network's radius approaches infinity, the system performance will, of course, degrade owing to the severe impact of large-scale path-loss. To mitigate such an ultra-long transmission, multiple gateways are one of the most effective solutions since they dramatically reduce the transmission distance (Tu et al., 2022b). Nonetheless, we leave such an interesting network for future work.

2 Other performance metrics such as average lifetime, and spectral efficiency can be effortlessly extended from the coverage probability (Nguyen et al., 2023b). Nonetheless, to comprehensively study these metrics, we leave it for the future work.

3 It is noted that the derived framework can be considered as an upper bound of the exact framework (Tu et al., 2020).

Appendix A.

Proof of Lemma 3.1

Let us begin the proof with the definition of the average transmit power as follows: (A1) $\begin{array}{rl}\mathbb{E}\left\{P_{i,k}\right\}& =\sum _{i=1}^N\frac{\chi P_P{\left|h_{i,o}\right|}^2{T}_{\mathrm{sleep},k}}{max\left\{1,{\Vert {\mathbf{p}}_o-{\mathbf{p}}_i\Vert }^{\alpha }\right\}{T}_{\mathrm{trans},k}}\le {\mathcal{P}}_k=\sum _{i=1}^N\frac{\chi P_P{\left|h_{i,o}\right|}^2{T}_{\mathrm{sleep},k}}{max\left\{1,t_o^{\alpha }\right\}{T}_{\mathrm{trans},k}}\\ & =N\frac{\chi P_P{T}_{\mathrm{sleep},k}}{T_{\mathrm{trans},k}}\mathbb{E}\left\{X=\frac{{\left|h_{i,o}\right|}^2}{max\left\{1,t_o^{\alpha }\right\}}\right\}.\end{array}$(A1) Here ${\mathbf{p}}_o$ and ${\mathbf{p}}_i$ are the position of the kth ED and ith power beacon, $\Vert .\Vert$ is the Euclidean distance and random variable (RV) $t_o$ refers to $t_o=|r_o-R/2|$. To obtain the mean of a random variable X, one first needs to obtain its cumulative distribution function (CDF) and probability density function (PDF) denoted by $F_X\left(x\right)$ and $f_X\left(x\right)$ that are given as follows: (A2) $\begin{array}{rl}{F}_X\left(x\right)& =Pr\left\{X=\frac{{\left|h_{i,o}\right|}^2}{max\left\{1,t_o^{\alpha }\right\}}\le x\right\}\stackrel{\left(a\right)}{=}{\int }_0^{R/2}\left(1-\exp \left(-\frac{xmax\left\{1,t^{\alpha }\right\}}{\mathrm{\Omega }}\right)\right)\frac{2}{R}\mathrm{d}t\\ & =1-\frac{2}{R}{\int }_0^{R/2}\exp \left(-\frac{xmax\left\{1,t^{\alpha }\right\}}{\mathrm{\Omega }}\right)\mathrm{d}t\\ & =1-\frac{2}{R}\left[{\int }_0^1\exp \left(-\frac{x}{\mathrm{\Omega }}\right)\mathrm{d}t+{\int }_1^{R/2}\exp \left(-\frac{x{t}^{\alpha }}{\mathrm{\Omega }}\right)\mathrm{d}t\right]\\ & \stackrel{\left(b\right)}{=}1-\frac{2}{R}\exp \left(-\frac{x}{\mathrm{\Omega }}\right)-\frac{2}{R}\frac{1}{\alpha }{\left(\frac{\mathrm{\Omega }}{x}\right)}^{\frac{1}{\alpha }}\left[\mathrm{\Gamma }\left(\frac{1}{\alpha },\frac{x}{\mathrm{\Omega }}\right)-\mathrm{\Gamma }\left(\frac{1}{\alpha },\frac{x}{\mathrm{\Omega }}{\left(\frac{R}{2}\right)}^{\alpha }\right)\right],\end{array}$(A2) where $\mathrm{\Gamma }\left(.,.\right)$ is the upper incomplete Gamma function (Gradshteyn & Ryzhik, 2007, (8.350.2)); $\left(a\right)$ is held by applying the CDF and PDF of the ${\left|h_{i,o}\right|}^2$ while $\left(b\right)$ is attained by computing the integration. The PDF of X is then derived by taking the first-order derivative with respect to (w.r.t.) to x as follows: (A3) $\begin{array}{rl}{f}_X\left(x\right)& =\frac{\mathrm{d}{F}_X\left(x\right)}{\mathrm{d}x}=\frac{2}{R\mathrm{\Omega }}\exp \left(-\frac{x}{\mathrm{\Omega }}\right)+\frac{2}{R}{\left(\frac{1}{\alpha }\right)}^2{\mathrm{\Omega }}^{\frac{1}{\alpha }}{\left(x\right)}^{-\frac{1}{\alpha }-1}\\ & \times \left[\mathrm{\Gamma }\left(\frac{1}{\alpha },\frac{x}{\mathrm{\Omega }}\right)-\mathrm{\Gamma }\left(\frac{1}{\alpha },\frac{x}{\mathrm{\Omega }}{\left(\frac{R}{2}\right)}^{\alpha }\right)\right]\\ & +\frac{2}{R}\frac{1}{\mathrm{\alpha x}}\left[\exp \left(-\frac{x}{\mathrm{\Omega }}\right)-\left(\frac{R}{2}\right)\exp \left(-\frac{x}{\mathrm{\Omega }}{\left(\frac{R}{2}\right)}^{\alpha }\right)\right].\end{array}$(A3) The mean of X is then computed as follows: (A4) $\begin{array}{rl}\mathbb{E}\left\{X\right\}& ={\int }_0^{\mathrm{\infty }}x{f}_X\left(x\right)\mathrm{d}x={\int }_0^{\mathrm{\infty }}x[\frac{2}{R\mathrm{\Omega }}\exp \left(-\frac{x}{\mathrm{\Omega }}\right)+x^{-1-\frac{1}{\alpha }}\frac{2}{R}{\left(\frac{1}{\alpha }\right)}^2{\left(\mathrm{\Omega }\right)}^{\frac{1}{\alpha }}\\ & \times \left[\mathrm{\Gamma }\left(\frac{1}{\alpha },\frac{x}{\mathrm{\Omega }}\right)-\mathrm{\Gamma }\left(\frac{1}{\alpha },\frac{x}{\mathrm{\Omega }}{\left(\frac{R}{2}\right)}^{\alpha }\right)\right]+\frac{2}{R}\frac{1}{\mathrm{\alpha x}}\\ & \times \left[\exp \left(-\frac{x}{\mathrm{\Omega }}\right)-\left(\frac{R}{2}\right)\exp \left(-\frac{x}{\mathrm{\Omega }}{\left(\frac{R}{2}\right)}^{\alpha }\right)\right]]\mathrm{d}x\\ & =\frac{2\mathrm{\Omega }}{R}\left(1+\frac{1}{\alpha -1}\left(1-{\left(\frac{R}{2}\right)}^{1-a}\right)\right).\end{array}$(A4) Here the integration related to the incomplete Gamma function is attained with the help of (Mathematica, n.d.) and we close the proof here.

Appendix B.

Proof of Lemma 3.2

Let us begin the proof with the definition of the CCDF as follows: (A5) $\begin{array}{rl}{\overline{F}}_Y\left(y\right)& =Pr\{\frac{{\left|h\right|}^2}{max\left\{1,d^{\alpha }\right\}}\ge y\}\stackrel{\left(a\right)}{=}{\int }_0^R\exp \left(-\frac{ymax\left\{1,d^{\alpha }\right\}}{\mathrm{\Omega }}\right)\frac{2d}{R^2}\mathrm{d}d\\ & ={\int }_0^1\exp \left(-\frac{y}{\mathrm{\Omega }}\right)\frac{2d}{R^2}\mathrm{d}d+{\int }_1^R\exp \left(-\frac{y{d}^{\alpha }}{\mathrm{\Omega }}\right)\frac{2d}{R^2}\mathrm{d}d\\ & \stackrel{\left(b\right)}{=}\frac{1}{R^2}\exp \left(-\frac{y}{\mathrm{\Omega }}\right)+\frac{1}{\alpha }\frac{2}{R^2}{\int }_1^{R^{\alpha }}{t}^{\frac{2}{\alpha }-1}\exp \left(-\frac{\mathrm{yt}}{\mathrm{\Omega }}\right)\mathrm{d}t\\ & =\frac{1}{R^2}\exp \left(-\frac{y}{\mathrm{\Omega }}\right)+\frac{1}{\alpha }\frac{2}{R^2}{\left(\frac{y}{\mathrm{\Omega }}\right)}^{-\frac{2}{\alpha }}\left[\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{y}{\mathrm{\Omega }}\right)-\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{y{R}^{\alpha }}{\mathrm{\Omega }}\right)\right].\end{array}$(A5) Here $\left(a\right)$ is attained by employing the CCDF of the exponential distribution and the PDF of a random point inside the disk while $\left(b\right)$ is held by changing variable $t=d^{\alpha }$ and the last equation is achieved by using the definition of the upper incomplete Gamma function. The CDF is defined as the complement of the CCDF and is evaluated as (A6) $\begin{array}{rl}{F}_Y\left(y\right)& =1-{\overline{F}}_Y\left(y\right)=1-\frac{1}{R^2}\exp \left(-\frac{y}{\mathrm{\Omega }}\right)\\ & -\frac{1}{\alpha }\frac{2}{R^2}{\left(\frac{y}{\mathrm{\Omega }}\right)}^{-\frac{2}{\alpha }}\left[\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{y}{\mathrm{\Omega }}\right)-\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{y{R}^{\alpha }}{\mathrm{\Omega }}\right)\right].\end{array}$(A6) The PDF is obtained by taking the first-order derivative w.r.t y as follows: (A7) $\begin{array}{rl}{f}_Y\left(y\right)& =\frac{\mathrm{d}{F}_Y\left(y\right)}{\mathrm{d}x}=\frac{1}{\mathrm{\Omega }{R}^2}\exp \left(-\frac{y}{\mathrm{\Omega }}\right)+\frac{4}{{\alpha }^2{R}^2}{\left(\frac{1}{\mathrm{\Omega }}\right)}^{-\frac{2}{\alpha }}\left(y^{-1-\frac{2}{\alpha }}\right)\\ & \times \left[\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{y}{\mathrm{\Omega }}\right)-\mathrm{\Gamma }\left(\frac{2}{\alpha },\frac{y{R}^{\alpha }}{\mathrm{\Omega }}\right)\right]+\frac{2}{\alpha R^2}\frac{1}{y}\\ & \times \left[\exp \left(-\frac{y}{\mathrm{\Omega }}\right)-R^2\exp \left(-\frac{y{R}^{\alpha }}{\mathrm{\Omega }}\right)\right].\end{array}$(A7) We complete the proof here.