Performance Analysis and Optimization of the Coverage Probability in Dual Hop LoRa Networks With Different Fading Channels

In this work, the performance evaluation and the optimization of dual-hop LoRa network are investigated. In particular, the coverage probability (Pcov) of edge end-devices (EDs) is computed in closed-form expressions under various fading channels, i.e., Nakagami- $m$ and Rayleigh fading. The Pcov under Nakagami- $m$ fading is computed in the approximated closed-form expressions; the Pcov under Rayleigh fading, on the other hand, is calculated in the exact closed-form expressions. In addition, we also investigate the impact of different kinds of interference on the performance of the Pcov, i.e., intra-SF interference, inter-SF interference (or capture effect) and both intra- and inter-SF interference. Our findings show that the impact of imperfect orthogonality is not non-negligible, along with the intra-SF interference. Moreover, based on the proposed mathematical framework, we formulate an optimization problem, which finds the optimal location of the relay to maximize the coverage probability. Since it is a mixed integer program with a non-convex objective function, we decompose the original problem with discrete optimization variables into sub-problems with a convex feasible set. After that, each sub-problem is effectively solved by utilizing the gradient descent approach. Monte Carlo simulations are supplied to verify the correctness of our mathematical framework. In addition, the results manifest that our proposed optimization algorithm converges rapidly, and the coverage probability is significantly improved when the location of relay is optimized.


I. INTRODUCTION
It is expected that there will be over 50 billion devices connecting to the Internet by the end of the year [1], making the Internet-of-Things (IoTs) a major component in the telecommunications industry. To support such a massive number of end-devices (EDs) networks, a couple of available technologies on the telecommunications market are taken into consideration. As the first, cellular networks are believed to The associate editor coordinating the review of this manuscript and approving it for publication was Yilun Shang . represent a suitable candidate owing to ultra-dense deployment of the base stations (BSs) and the well-established standards. With ultra-dense deployment, however, the BSs' power consumption has become one of the most significant issues in the information and communications technology (ICT) field, which accounts for approximately 2% of worldwide CO 2 emissions. In addition, with high capital expenditure (CAPEX) and operational expenditure (OPEX), the mobile networks seem not to be a wise choice for this kind of network. As a result, low power wide area networks (LPWAN) are regarded as the suitable technology for the VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ massive IoTs networks. The main advantage of LPWAN is that it is able to connect enormous low-power end-devices with a simple protocol and infrastructure [2]. Among all available LPWAN technologies, i.e., SigFox, Ingenu and long-range (LoRa), LoRa is emerging as the most promising and attracts the attention of numerous researchers from both industry and academia. The main reason for the success of LoRa is that LoRa's signals are modulated by using chirp spread spectrum (CSS) modulation instead of the conventional modulations, i.e., QAM, PSK and FSK, which is proven to better resist fading and noise. Furthermore, a wide set of parameters also contributes to the popularity of LoRa among all LPWAN technologies. In particular, by actively adjusting the spreading factor (SF), the transmit power, and the bandwidth (BW), LoRa is capable of satisfying long-range transmissions with different requirements and low power consumption. Moreover, to maximize the performance of the entire networks, LoRa separates its coverage area into a set of non-overlapping regions, where each region will be assigned a unique value of SF and transmit power: in particular, the nearer the EDs is, the lower the SF and the smaller the transmit power. The aim of this resource allocation is not only to reduce the inter-SF interference but also to conserve power consumption by the EDs. Although smart resource allocation is yielded, the edge EDs still suffer from significant interference compared with end-devices around the gateway. This problem is even more serious in LoRa due to the lack of power control at uplink transmission, or any kinds of reliable signaling protocols. Consequently, the performance of edge-ED apparently becomes the bottleneck of the entire networks.
On the other hand, relaying communications is proposed and regarded as an efficient way to improve the performance in wireless networks [3], [4]. It is evident that with the help of relay, the transmission distance is shortened: hence, the reliability, of course, is ameliorated and the transmit power is dramatically declined. Moreover, another benefit of relaying communications is extension of the coverage area of the networks. Thus, in this work, we study and enhance the performance of edge-ED in LoRa networks with the aid of relaying. Before briefly summarizing the main contributions and novelties of the present paper, some state-of-the-art LoRa networks with and without the help of relay are first reported in the sequel.
The performance evaluation of LoRa networks was studied under different circumstances. In [5], the coverage probability (Pcov) was studied under the assumption that inter-SF interference was absent. However, as pointed out in [6]- [8], the inter-SF interference or called capture effect could not be ignored in practical scenarios. In [8], the Pcov, which was considered with respect to both the intra-and inter-SF interference, as well as interference from different technologies, was investigated. Nevertheless, the metric was computed via numerical computation, or no closed-form expression was provided. In [6], antenna diversity was used to enhance the performance of LoRa networks. However, multiple antennas at either gateway and/or EDs seems to be impractical in LoRa networks owing to low cost transceiver. The ergodic capacity was studied in [9].
The performance of relay networks were studied in [10]- [12]. Particularly, the outage performance of dual-hop with amplify-and-forward (AF) fixed gain relaying was studied in [10]. The ergodic capacity of multi-hop decodeand-forward (DF) relaying was studied in [12]. The results showed that the optimal rate adaption attained the highest spectral efficiency. The combination of cooperative networks with other techniques, i.e., simultaneous wireless information and power transfer (SWIPT), non-orthogonal multiple access (NOMA) were investigated in [13]- [15]. In [13] the throughput of the cooperative networks where the relay operated based on the harvested energy was investigated by using AF protocol. The outage probability of the SWIPT-enabled cooperative networks where the positions of the relays were randomly located was studied in [14]. The results illustrated that under some assumptions, the SWIPT-enabled cooperative networks was able to achieve the same diversity gain as the conventional cooperative networks. The throughput and outage probability of the cooperative SWIPT NOMA networks were addressed in [15]. The results demonstrated that the use of SWIPT did not jeopardize the diversity gain compared to the conventional NOMA. In addition, relaying communications was also widely utilized in different networks/topics, such as cognitive radio networks [16], [17], physical layer security, device-to-device (D2D) and cellular networks [20], [21]. The secrecy performance with relay selection under impact of co-channel interference was investigated in [18] and the maximum capacity of relay-aided D2D communications was studied in [19]. However, the application of relaying communications into LoRa networks remains in its infancy. In [22], a multi-hop concurrent transmission LoRa network was investigated. The main target of this work was to study the impact of concurrent transmission on the performance at the LoRa receiver, i.e., capture effects, energy spreading effect and so forth. Reference [23] studied the performance of the practical dual-hop LoRa networks. Particularly, the experiment was deployed at the campus of the university of Pau, France where the battery-based relay was placed between end-devices and gateway in order to help in forwarding end-devices' packet to gateway. The paper, however, mainly focused on the power consumption at the relay by calculating the power consumption and then proposing a wake-up algorithm in order to save power consumption.
In this work, in contrast with these above-mentioned works, we focus on the performance of edge-EDs in dual-hop decode-and-forward relaying in LoRa networks under different types of fading channels. In particular, in contrast with [22], we take into consideration the scenario where the transmission is divided into different time-slots to maintain the interference to be as minimal as possible and to improve the performance of edge-ED. Compared with [23], we are interested in studying the performance of ED instead of the relay. The main contributions and novelties are summarized as follows: • The Pcov of edge-ED is computed in closed-form expressions under different fading channels, i.e., Nakagami-m and Rayleigh. To be more specific, the Pcov under Nakagami-m is computed in the approximated but closed-form expression, while the framework under Rayleigh fading is calculated in the exact closed-form expression.
• We address all types of interference in LoRa networks, i.e., intra-SF, inter-SF and both intra-and inter-SF interference. Numerical results show that the impact of capture effects on the performance of LoRa networks is not non-negligible, along with the intra-SF interference.
• An effective algorithm is provided to optimize the position of the relay, which maximizes the Pcov of the edge-EDs.
• The baseline system without relay is also presented to highlight the benefits of the proposed networks. The rest of this paper is organized as follows: In Section II, the system model is introduced. The performance of the Pcov is computed in Section III. Also in this section, an optimization problem of the Pcov with respect to the position of the relay is formulated and solved by an effective algorithm. In Section IV, numerical results based on Monte Carlo simulations are provided to confirm the correctness of the proposed frameworks. Finally, Section V concludes this paper.
Notations: Main notations and mathematical symbols/shorthand are provided in Table 1.

II. SYSTEM MODEL
Let us consider uplink LoRa networks where the desired end-device denoted by S communicates to the gateway, G, via a relay denoted by R. We assume that S always has packets to send to the gateway. In LoRa, depending on the distance from the ED to the gateway, an appropriate spreading factor (SF) as well as the transmit power of EDs will be assigned in order to guarantee fairness among EDs at different locations. To be more specific, the considered network is split into six non-overlapping regions with equal distance denoted by SFk, k ∈ {7, . . . , 12}, and both the spreading factor and transmit power are assigned based on the incremental rule: the closer the gateway, the smaller the SF and transmit power, as shown in Fig. 1.
In addition to S, R and G, the considered networks also comprise N = 12 k=7 N k EDs, which act as interferer to the intended link; N k is the number of interferers from region SFk. Interference from other technologies which operate at the same industrial, scientific and medical (ISM) band is not considered [5]. The position of the relay, denoted by V R = v x , v y , is assumed to locate at SFo, o ∈ {7, . . . , 12} and is changeable, while the locations of both S and G are assumed to be fixed. Here, v x , v y are the horizontal and vertical coordinates of relay R, respectively. For simplicity, we assume that G is fixed at the origin and that S is at distance R from G, located in the furthest region, i.e., SF12, as shown in Fig. 1.
The transmission from S to G occurs in two consecutive time-slots or two phases. In the first phase, the ED of interest transmits its signals to the relay, and the signals received at the relay are formulated as: where h S,R , h k i,R are the channel coefficient from S and interferer i of region SFk to relay R, respectively, and follow the Nakagami-m distribution with shape and spread parameters denoted by m S,R , m k i,R θ S,R and θ k i,R , correspondingly. The channel gain, denoted by h S,R 2 , h k i,R 2 , follows a Gamma distribution with the following shape and scale parameters m S,R , m k i,R and 1/β S,R , are the corresponding large-scale fading from S and interferer i of SFk to R including shadowing; d X ,Y is the Euclidean distance from X to Y and is computed as Y } are the horizontal and vertical coordinates of node Z . In this work, time is slotted, and we further assume that the fading remains constant during one time-slot and changes between timeslots. P k is the transmit power of ED belonging to the region with SFk, and we assume that all EDs have the same transmit power in each region. x S , x k i,R are modulated signals of S and interferer i of SFk modulated by the patented CSS modulation with unit power, i.e., E |x is the activation function of interferer i of SFk and follows the Bernoulli distribution with the success probability is the bit rate of ED in region SFk (in bit/s), which is provided in Table 2. L pac , T in , CR and BW are the packet length (in bits), the interarrival time between two packets (in seconds), the coding rate and the transmission bandwidth (in Hz), respectively. In the present work, we assume that all EDs, regardless of region, have the same packet length, inter-arrival time, coding and transmission bandwidth. In (1), n R is the AWGN noise at the relay node with zero mean and variance [5]: where the first part is thermal noise normalized to 1 Hz. The second part, NF, is noise figure of the receiver (in dBm) and the last part contains the effects of the used bandwidth.
We notice that, in (1), the term represents the aggregate interference from the signals using the same SF or intra-SF interference and signals from different SFs or inter-SF interference.
At the end of the first phase, the relay R decodes the signal from S following re-modulating and forwarding the information to gateway G. In the present paper, the decode-andforward protocol is utilized. The core reason of this utilization is that compared to amplify-and-forward protocol, the DF protocol achieves better performance as well as requires less complexity hardware at the relay [24]. As a consequence, the received signals at the gateway, denoted by y G , are formulated as: where P R = P o , o ∈ {7, . . . , 12} is the transmit power of the relay; the explicit values of P o are available in Table 2; h R,G , h k i,G are the channel coefficients from R and interferer i of SFk to the gateway; x R , x k i,G are re-modulated signals of S at R and signals from interferer i of SFk; n G indicates AWGN noise at the gateway; χ k i,G is the activation function. In this paper, we assume that the active interference of all regions is exactly the same for two phases. The asymmetric case where the active EDs of the first and second phase are independent can be derived in straightforward fashion by using our following mathematical frameworks, since the impacts of interferer EDs on relay and gateway are noncorrelated.
Under the considered networks, the signal-to-noise ratio (SNR) of the transmitted signals from X to Y , denoted by SNR XY , is formulated as follows where P X is the transmit power of node X ; σ 2 Y is the noise variance at receiver Y and h X ,Y 2 is the channel gain from X to Y . The signal-to-interference ratio (SIR) of the packets sent from node X of region SFo, o ∈ {7, . . . , 12} to node Y impaired by interference from SFk, k ∈ {7, . . . , 12}, is formulated as follows where, N k = p k A N k is the number of active EDs belonging to SFk; . is the ceiling function; h k i,Y 2 , P k are the channel gain and transmit power of interferer i from SFk to receiver Y .
In the sequel, the coverage probability of the transmission from S to G is computed in closed-form expressions under different scenarios.

III. COVERAGE PROBABILITY ANALYSIS
In LoRa, the coverage probability refers to the probability that an arbitrary ED is in coverage or that its packets are successfully transmitted to the gateway. To be more specific, one packet which operated at SFo, o ∈ {7, . . . , 12}, is considered to be decoded correctly if the two following conditions are satisfied simultaneously: i) its SNR is greater than a given threshold, q o , where q o values are provided in Table 2; ii) its SIR versus other packets from the same or different SFs k are larger than the rejection threshold, o,k (in dB), o, k ∈ {7, . . . , 12}, [5], [6]. Here, o,k is the o-row and k-column entry of matrix , which is provided as follows [26, For example, if one packet is sent at SF9, then it can be decoded error-free, provided that its SIR versus the packet from SF7 is not less than −15 dB and that its SIR versus the packet from SF9 is at least 1 dB.
To simplify the mathematical formulas we set ϑ S,R = q S,R β S,R σ 2 R P S . In the sequel, two following useful Lemmas are proposed in order to compute the coverage probability and are given as where (.), γ (., .) are the gamma function and lower incomplete gamma function; ς and ξ are provided as follows Proof: The proof is available in Appendix . Lemma 2: Given two Gamma RVs X and Y with corresponding shape and scale parameters, The proof is available in Appendix . Based on the outcomes of Lemmas 1 and 2, the coverage probability (Pcov) of three different interference circumstances, namely, intra-SF interference, inter-SF interference and both intra-SF and inter-SF interference, are formulated and computed by three following Theorems.
Theorem 1: Assuming that relay R is located at region SFo, o ∈ {7, . . . , 12}, under intra-SF interference the coverage probability of the signals from S to R and from R to G is formulated in (10): and it is computed in (11), as shown at the bottom of the next page. Here C 1 (q o ) is the probability that the SNR values of both hops are greater than a given threshold. It is noted that since node S is always located at the edge of the networks or SF12, the SNR threshold, q S,R , is subsequently always q 12 , q S,R = q 12 ; q o , on the contrary, is changeable because of the flexible position of the relay. C intra 2 is the probability that the SIR values of two hops are greater than the rejection threshold, intra X,Y , X ∈ {S, R}, Y ∈ {R, G}, and rely on the region where the packet is sent; in addition, intra X,Y is also the diagonal element of matrix in (6). For example, the rejection threshold of packets sent by S, which is located VOLUME 8, 2020 at SF12, is intra S,R = 12,12 ; the same for the second hop, we have intra R,G = o,o . Proof: The proof is available in Appendix . Theorem 2: When the inter-SF interference is taken into consideration, the coverage probability of S with the help of relay R located at region SFo, o ∈ {7, . . . , 12} is formulated in (12): and it is computed in (13), as shown at the bottom of this page. Here and computed in (15), as shown at the bottom of this page.
Here, C 1 q y is the same as Theorems 1 and 2. C both 2 is defined as the probability that both SIR values under both inter-and intra-SF interference from S to R and from R to G are larger than the rejection threshold. In this context, the SF values of the intended packet and interference packet are not necessarily the same. In particular, both S,R = 12,k , k ∈ {7, . . . , 12} and both R,G = o,k ; o, k ∈ {7, . . . , 12}. Proof: The proof can be intuitively derived by combining the findings from Theorems 1 and 2. Particularly, C 1 (q o ) is obtained from (47), and C both 2 comes from the multiplication of (49) with (52). We close the proof here.
Remark 1: By directly inspecting (15), it is apparent that although the Pcov can be computed in closed-form expression. The results, however, are solely approximated due to the approximation of the aggregate interference. Thus, Corollaries 1 and 2 are provided not only to simplify the mathematical framework but also to obtain the exact closed-form expression where Rayleigh fading is taken into consideration.
To be specific, the coverage probability under Rayleigh fading of all scenarios, i.e., intra-SF, inter-SF and all interference, are computed in the exact closed-form expressions. The derivation, nevertheless, provides solely the case of all interference, while the two remaining case studies are directly attained from the general case.
and it is computed in (17), as shown at the bottom of the next page.
Here ω X ,Y , X ∈ {i, S, R}, Y ∈ {R, G}, in (17) is derived from β X ,Y by letting m X ,Y = 1, and we have ω X ,Y = L X ,Y θ X ,Y . Proof: The proof is available in Appendix . Corollary 2: The Pcov of intra-SF and inter-SF interference under Rayleigh fading are provided in (18) and (19), as shown at the bottom of the next page.
In the next section, we will maximize the performance of the coverage probability by optimizing the relay position under the considered networks. To be more specific, the cov- erage probability of all intra-and inter-SF interference under Rayleigh fading are yielded.

A. COVERAGE PROBABILITY MAXIMIZATION
We now formulate an optimization problem that maximizes the coverage probability in (17) with respect to the location of the relay in the different zones as 1 Since problem (20) is non-convex and the objective function should be redefined as a function of each zone, a standard convex optimization toolbox is not applicable for solving this problem. For further manipulation, we recast problem (20) for a given zone o as where P o v x , v y is a coverage probability for the particular zone o, which is derived from the general formula (17) and presented in the closed-form expression in (25) with (20) is formulated based on the optimistic assumptions, for example the perfect knowledge of interferers and is only solved offline. However, the solution to this problem can be used together with deep neural network for online resource allocation, similar to what has been done in [27], and we leave this for the future work.
Even though problem (24) is still non-convex, its feasible domain is convex and the objective function is continuous. Consequently, according to Weierstrass' theorem, a globally optimal solution exists [28], [29]. In order to find a local optimal solution, we now introduce the partial Lagrangian function of problem (24) as which is based on the fact that maximizing P o v x , v y is equivalent to minimizing −P o v x , v y . Despite the inherent non-convexity, we can take the first-order derivative of the partial Lagrangian function with respect to v x and v y as ∂L and the closed-form expressions are shown in (29) and (30), as shown at the bottom of the next page, with the other supported variables defined as: Here, we notice thatḟ (z) = ∂f /∂z,ġ (z) = ∂g/∂z is the partial derivative of f and g with respect to z, z ∈ v x , v y . From initial values v (0) x and v (0) y in the feasible domain, and by exploiting gradient descent to find a local minimum [30], the coordinate is first updated at iteration n aŝ where the step size τ > 0 is sufficiently large for the direction of steepest descent. Due to the constraints on the zone in (24), the coordinate of the relay should be updated by checking the boundary values as After a number of iterations, the updates in (34) and (35) converge to a local solution for zone o. The stopping criterion can be defined based on the variety between two consecutive iterations, which should be less than a given value: By assuming that the convergence holds at iteration n, we obtain a local solution as for which the locally optimal objective value is P . We stress that the above optimization procedure is applied to a particular zone, and it will be repeated for all remaining zones to find a good relay position. By gathering all of the optimized solutions of all six considered zones, the solution to problem (20) is obtained by computing: Our proposed optimization approach to maximize the coverage probability is summarized in Algorithm 1. 2

B. BASELINE PERFORMANCE
In this section, the performance of LoRa networks without the help of relay is provided. Particularly, the coverage probability (Pcov) of direct transmission under both Nakagami-m and Rayleigh fading are provided by Proposition 1. It is noted that 2 Algorithm 1 obtains a local optimum to problem (24) with low computational complexity in each iteration thanks to the use of gradient descent, which only needs to evaluate the first derivative of the partial Lagrangian function with respect to each optimization variable. The performance of Algorithm 1 compared to the other benchmarks are presented by numerical results in Section IV.

Algorithm 1 A Local Solution to Problem (20) by Using Gradient Descent
Input: K 0 , (x S , y S ), (x G , y G ), η, θ S,R , θ R,G , P S , q S,R , σ 2 S , σ 2 R , q o , P R , 12,k (27) and (28) . Then, obtain the optimized coordinates of the relay by solving (38). Output: The coordinates of the relay v * x and v * y .
although the coverage probability of LoRa networks under Rayleigh fading was well-studied in [5], the performance under Nakagami-m is still missing in the literature. As a result, the main purpose of this baseline system is to not only act as an benchmark against of proposed relay-aided LoRa networks but also highlight the impact of the fading parameter, m, on the performance of edge-node LoRa networks. The coverage probability (Pcov) of direct transmission is provided by the following Proposition.
Proof: The proof is derived easily by following the same stages in Theorem 3.

IV. NUMERICAL RESULTS
In this section, numerical results are provided to verify the correctness of our mathematical frameworks. Particularly, the following setups are used: BW = 250 KHz, NF = 6 dBm, η = 3, f c = 868 MHz, L pac = 10 bytes, CR = 4/5, T in = 60s; we consider the rectangular networks with radius R = 6000 m; the horizontal and vertical networks are from 0 → R and −R/2 → R/2, respectively; and v x , v y = (2.5R/6, −R/12), (c S , w S ) = (R, 0) and (c G , w G ) = (0, 0). The transmit power of interference EDs, as well as relay, are changeable and rely on the located region, while the transmit power of source node, P S , is always fixed and equal to P 12 = 17 dBm; the detailed values of the transmit power of all regions are available in Table 2. The SNR threshold, q o , o ∈ {7, . . . , 12}, is also given in Table 2. As for the small-scale fading, the following parameters are considered: Then, the set of active interferer of each region is identified based on the Bernoulli distribution with its corresponding success probability p k A , k ∈ {7, . . . , 12}. This set of active interferer remains active for the whole transmission which comprises of two phases. In the first phase, the source node transmit packets to the relay, the relay decodes, then forwards to the gateway at the second phase. The performance metric is computed for each packet based on the SNRs and SIRs at the received nodes, i.e., relay and gateway. For each realization, 200 packets are transmitted from the source node and the procedure is iterated for 10000 realizations of different locations of the interferer. Thus, the Pcov is averaged over 2 million packets. Fig. 2 shows the coverage probability with respect to the SNR threshold q S,R = q S,G with all kinds of interference, i.e., intra-SF, inter-SF and both interferences under Nakagami-m distribution. It is obvious that our mathematical frameworks exactly overlap with respect to Monte Carlo simulations, thus verifying the correctness of our derivation. Firstly, we observe that the larger the q S,R = q S,G , the smaller the Pcov. This can be justified in straightforward fashion by direct inspection of the definition of the Pcov. Moreover, there is no doubt that the performance of the coverage probability of the edge-node such as S is significantly ameliorated with the help of the relay, e.g., approximately 0.1 when q S,R = q S,G = −20 dBm. It emphasizes that under some specific value of q S,R = q S,G , i.e., q S,R = q S,G = −4 dBm, the proposed systems can improve the Pcov of edge-node up to 90% compared with the conventional single-hop LoRa networks. It is also evident that with the same value of the Pcov, the proposed dual-hop LoRa networks can support higher SNR threshold than the baseline systems, e.g., approximately 10 dBm when Pcov = 0.6. It means that we are able to enhance the QoS at the edge-node without degrading the coverage probability.
In addition, Fig. 2 also confirms the necessity of taking into consideration the impact of the capture effect or imperfect orthogonality in LoRa networks. Particularly, under the degradation of both intra and inter-SF interference, the Pcov slightly decreases when q S,R = q S,G is relatively small. If q S,R = q S,G continues to increase, then all curves become indistinguishable due to the change of networks from interference-limited to noise-limited scenario. This can be explicated by directly inspection (15), when q is small, the Pcov is dominated by the term 12 k=7 J 1 12,k β S,R P 12 or the system is under interference-limited and when q is sufficiently large, Pcov changes to noise-limited and dominated by . Hence, there is no differences between intra-, inter and both interferences in this region. Fig. 3 illustrates the behavior of Pcov with respect to the number of interferer, N , and different small-scale fading, i.e., Rayleigh, Nakagami-m and no fading, respectively. Observing this figure, we see that when number of interferer is relatively small or the system is under sparsely-loaded scenario, the no fading, of course, is the best followed by Nakagami-m and Rayleigh fading is the worst. This can be explained that when N is small, the system is in noise-limited regime or the impact of the aggregate interference is minority compared to the AWGN noise. Under this context, the no fading outperforms others. On the other hand, when the network is densified or in fully-loaded scenario, N 1,  20]; for Nakagami-m, we use the typical setup except for θ S,R = θ R,G = 2; θ S,G = 8 in order to have equal channel gain of all scenarios. Solid lines are plotted by using (17) and (15). Markers are from Monte-Carlo simulations.
the network is under interference-limited regime and the no fading becomes the worst and Rayleigh fading is the best one. This phenomenon can be easily explained that when keep increasing number of interferer N , the aggregate interference of case no fading will increase faster than Nakagami-m fading and Rayleigh due to its better channel gain to the gateway. The intended link, on the contrary, is constant with this augmentation of the interference. As a result, the SIRs of the no fading declines fastest among three case studies followed by the Nakagami-m fading, thus, the Pcov of both the no fading and the Nakagami-m decrease quicker and become worse than the Rayleigh fading. It is obvious that network densification is monotonic decreasing the performance of the LoRa networks. However, if the considered metric is the potential area spectral efficiency (PSE) [31], we believe that there exists an optimal value of N that maximize PSE, studying the PSE, nonetheless, is out of scope of the current paper and is left to the future work. Again, this figure justifies the accuracy of our mathematical frameworks against results based on Monte Carlo simulation. Furthermore, it is interesting that the Pcov exhibits a down-stair property with respect to N . This trend can be obviously explicated: the number of active EDs of each region is a ceiling function rather than a true continuous function, i.e., N k = p k A N k . This means that if N k does not increase to a sufficiently large value, the N k , as a result, remains constant or changes with an epsilon pace.  (17) and (40). Markers are from Monte Carlo simulations. Fig. 4 shows the performances of coverage probability for both the proposed and baseline frameworks under Rayleigh fading with different values of the path-loss exponent. We observe that increasing the path-loss exponent, η, will, of course, decrease the coverage probability as the large-scale path-loss monotonically increases with path-loss exponent. Moreover, as pointed in Fig. 3, it is not surprising that the Pcov also follows the stair-shape when N keeps increasing. The results also consolidate the essential nature of the aid of the relay node, especially with respect to dense networks or N 1, as the gap between the ''Pro'' and the ''BL'' curves is approximately 0.3. Additionally, the curves with small path-loss exponent fall swiftly compared to large path-loss exponent. The principal reason for this trend is that under large path-loss exponent scenario, the interference  by edge nodes or EDs which are away from the receiver is negligible so the Pcov decreases with lower pace compared to case small path-loss exponent.
We show the CDF of the coverage probability by utilizing the different methods to locate the relay in Fig. 5. Locating the relay randomly in the network area results in the worst coverage probability of 0.43 on average. This is because bad locations for the relay may be encountered as a consequence of the relay's location randomness. Subsequently, by selecting a good heuristic fixed location for the relay (i.e., v x = 2500 and v y = −500, the coverage probability is 0.70. However, the fixed location is not always a good option when the receiver changes its location. Therefore, it calls for carefully relocating the relay's location whenever the receiver is relocated. Algorithm 1 yields the best coverage probability of 0.86 on average. It hence manifests superior improvements when optimizing the relay location. Fig. 6 plots the convergence of Algorithm 1 for all considered zones. The results numerically verify the convergence of our proposed solver, which requires less than 200 iterations to reach a fixed point for all the six zones. Compared with the initial point, optimizing the relay location in each zone achieves 7% to 30% improvement in the coverage probability. In addition, Fig. 6 visualizes that the convergence properties of each zone differ from each other. The convergence of regions 8 and 9 exhibits fluctuations along iterations, while the remaining zones exhibit a monotonic convergence property. Moreover, the coverage probability in zones 9 and 10 are better than the remaining since the distance from the relay to the source and receiver is relatively equal to each other without any drastic loss.

V. CONCLUSION
In this work, the coverage probability of the edge node is studied via the aid of one relay node which is randomly distributed throughout the networks. All results are presented in closed-form expressions, where the results under Rayleigh fading, in addition, are also computed in exact closed-form. Our findings show that it is obvious that the quality of the edge node is significantly ameliorated by using the relay. Furthermore, if the position of the relay is properly placed, we even further improve the network performance. In addition, the current work can be extended to several directions. First, it is no doubt that studying the performance of the considered networks under generalized fading channel is an interesting direction. In fact, some well-known generalized fading distributions like κ − µ, α − µ and η − µ are better described the small-signal variations rather than Nakagami-m fading [32], [33]. Next, another possible direction is to apply Intelligent Reflecting Surface (IRSs) and compare it performance with the traditional relay-aided LoRa networks [34], [35]. It is also possible to combine the advantages of both mathematical-based modelling with the power of deep learning to optimize the LoRa networks, i.e., the throughput of whole networks under more practical conditions [36], [37].

PROOF OF LEMMA 1
In this section, the statistics of the approximated RV of sum of N independent and non-identical distributed (i.n.i.d) Gamma RVs denoted by X ≈ X = N i=1 X i is proven. First, it should be noted that no exact closed-form expressions of both CDF and PDF of X are available in the literature. We, as a result, approximate X by X , which also follows the Gamma distribution based on the moment matching method (MMM) [38]. Let us denote ς and ξ , as the shape and scale factors of X . According to the MMM method, the mean and variance of two RVs are equivalent and are provided as follows where α i , β i are the scale and shape parameters of RV X i ; E {.} and Var {.} are the expectation and variance operator. From (41), we have the following: Next, by substituting ς in (42) to (41), we are able to compute ξ as follows Next, by substituting ξ in (43) to (42), we obtain the ς as follows Finally, by substituting ς and ξ from (43) and (44) into the formula of the CDF and PDF of Gamma RV, we obtain (7) and close the proof here.

PROOF OF LEMMA 2
The CCDF of RV Z , that is, the ratio of two Gamma RVs X and Y , is computed as follows. First, let us formulate the CCDF as follows where (a) is achieved by using the definition of the CCDF, i.e., F X (x) = 1 − F X (x) = Pr {X ≥ x}; (b) is replacing the formula of both CCDF and PDF of X and Y , respectively. Finally, (c) is attained by using the results from [39, Eq. 6.455].
Here a, b, c, z) is the Gaussian Hypergeometric function. We finish the proof here.

PROOF OF THEOREM 1
From (10), we observe that the Pcov constitutes two parts, C 1 (q o ) and C intra 2 . Thus, let us begin computing the C 1 (q o ) as follows where (a) is immediately attained by using the fact that the channel gains between the first and second hops are independently combined with the use of the CCDF of the Gamma distribution [40, Table 5-2]; q S,R = q 12 and q o are provided in Table 1. Now, let us investigate the second condition, C intra 2 , as follows where (a) is obtained by using the outcomes of Lemma 2 with respect to the CCDF of the ratio of two Gamma RVs, |hS,R| 2 I SF12 and |hR,G| 2 I SFo . Here, the shape and scale parameter of Gamma RV I SFl , l ∈ {o, 12}, is obtained thanks to Lemma 1 and provided as follows where l ∈ {y, 12}. We complete our proof here.

PROOF OF THEOREM 2
The coverage probability of the inter-SF interference scenario is calculated in this section. Let us re-write the definition for this case study as follows = Pr SIR inter SR ≥ inter S,R , SIR inter RG ≥ inter R,G . (51) By direct inspection (51), it is evident that C 1 (q o ) is identical to C 1 (q o ) in the case of intra-SF interference. As a consequence, we intuitively achieve the expression of C 1 (q o ) provided in (47) We conclude the proof here.

PROOF OF COROLLARY 1
Let us begin this section by re-writing the definition of Pcov in (16) The C Ra 1 (q o ) can be derived effortlessly from (47) by letting m S,R = m R,G = 1, as Here, (a) is obtained by using an identity [40,Eq. 8.352]; ω X ,Y , X ∈ {S, R}, Y ∈ {R, G} is derived from β X ,Y by letting m X ,Y = 1, and we have ω X ,Y = L X ,Y θ X ,Y . As for the C Ra,both 2 , it is computed as follows where (a) is derived as where (a) is attained by using the definition of the CCDF; (b) is derived by using the definition of the moment generating function (MGF) of exponential RV with mean µ as follows The remaining probability is freely derived by following the same procedures as those of (57). Finally, by combining the outcomes of (56) and (56), we obtain (17) and finish the proof here.