Aggregation and Resource Scheduling in Machine-type Communication Networks: A Stochastic Geometry Approach
Onel L. Alcaraz López Hirley Alves Pedro Nardelli Matti Latva-aho

	
	
		

					
Aggregation and Resource Scheduling in Machine-type Communication Networks: A Stochastic Geometry Approach
    
		Onel L. Alcaraz Lpez, Student Member, IEEE, 
			Hirley Alves, Member, IEEE,
			Pedro H. J. Nardelli, 
			Matti Latva-aho, Senior Member, IEEE
				
		
		Onel L. Alcaraz Lpez, Hirley Alves, Matti Latva-aho are with the Centre for Wireless Communications (CWC), University of Oulu, Finland.{onel.alcarazlopez,hirley.alves,pedro.nardelli,matti.latva-aho}@oulu.fi
		
		Pedro H. J. Nardelli is with Laboratory of Control Engineering and Digital Systems, Lappeenranta University of Technology, Finland. pedro.nardelli@lut.fi

		This work is partially supported by Academy of Finland (Aka) (Grants n.303532, n.307492), SRC/Aka (n. 292854), and by the Finnish Funding Agency for Technology and Innovation (Tekes), Bittium Wireless, Keysight Technologies Finland, Kyynel, MediaTek Wireless, Nokia Solutions and Networks.
	
    
===================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================





Data aggregation is a promising approach to enable massive machine-type communication (mMTC). This paper focuses on the aggregation phase where a massive number of machine-type devices (MTDs) transmit to aggregators. By using non-orthogonal multiple access (NOMA) principles, we allow several MTDs to share the same orthogonal channel in our proposed hybrid access scheme. We develop an analytical framework based on stochastic geometry to investigate the system performance in terms of average success probability and average number of simultaneously served MTDs, under imperfect successive interference cancellation (SIC) at the aggregators, for two scheduling schemes: random resource scheduling (RRS) and channel-aware resource scheduling (CRS).  We identify the power constraints on the MTDs sharing the same channel to attain a fair coexistence with purely orthogonal multiple access (OMA) setups. Then, power control coefficients are found, so that these MTDs perform with similar reliability. We show that under high access demand, the hybrid scheme with CRS outperforms the OMA setup by simultaneously serving more MTDs with reduced power consumption.




§ INTRODUCTION


Machine-type Communication (MTC) is typically cited as an integral use case for fifth generation (5G) cellular networks, where fully automatic data generation, exchange, processing and actuation among intelligent machines, are on the agenda. With the rapid penetration of embedded devices, MTC is becoming the dominant communication paradigm for a wide range of emerging smart services including healthcare, manufacturing, utilities, consumer goods and transportation. Specifically, massive MTC (mMTC), as the name suggests, is about massive  access by a large number of devices, that is, about providing wireless connectivity to enormous number of often low-complexity low-power machine-type devices (MTDs). Of course, the access is a major concern and a flouring research area. Different methods have been proposed and addressed to provide more efficient access, e.g., access class barring <cit.>, prioritized random access <cit.>, backoff adjustment scheme <cit.>, delay-estimation based random
access <cit.>, distributed queuing <cit.>. Another promising way to deal with the massive connection problem comes from the concept of data aggregation. Instead of being directly connected to the access network, e.g., by equipping MTDs with own Subscriber Identity Module (SIM) card to have cellular connectivity, the idea is that MTDs may organize themselves locally, creating MTC area networks and exploiting short-range technologies. These MTC area networks may then connect to the core networks through MTC gateways or data aggregators. This is a key solution strategy to collect, process, and communicate data in MTC use cases with static devices, especially if the locations of the devices are known, such as smart utility meters <cit.> or video surveillance cameras <cit.>. The traffic from MTDs is first transmitted to the designed data aggregator, while the aggregator then relays the collected packets to the core network <cit.>, e.g., the base station (BS), thus, reducing the congestion and the power consumption at the MTD side <cit.>.

When aggregating such huge number of MTDs, the density of the aggregators, although considerable smaller than the density of the MTDs, obviously will still be large, and the interference coming from the devices sharing the same resource could be significant. There is limited literature that  characterizes the interference in mMTC with data aggregation. The analysis in <cit.> shows how the wireless channel, transmit power, and random deployment of data collectors affect the signal-to-interference-plus-noise ratio (SINR) distribution in random access networks with randomly deployed sensor nodes, and for some special cases, a simple form of signal-to-noise ratio (SNR) or SINR distribution is found. An energy-efficient data aggregation scheme for a hierarchical MTC network is proposed in <cit.>, while authors develop a coverage probability-based optimal data aggregation scheme for MTDs to minimize the energy density of the network.
In <cit.>, a theoretical and numerical framework, which aims to assess, model and characterize the network energy consumption profile, is presented by exploiting the stochastic geometry tools. 
By incorporating some intelligence in the aggregator, the network performance improves, as shown in <cit.>  
for resource scheduling strategies. Among them, only <cit.> considers a more realistic scenario with a multi-cell network, hence the inter-cell interference, which is a critical issue, is taken into account. 
The authors introduce a tractable two-phase network model for mMTC, where MTDs first transmit to their serving aggregators (aggregation phase) and then, the aggregated data is delivered to BSs (relaying phase). Key metrics such as the MTD success probability, average number of successful MTDs and probability of successful channel utilization are investigated for two resource scheduling schemes RRS,  where aggregators randomly allocate the limited resources to MTDs, and CRS, where aggregators allocate resources to the MTDs having better channel conditions. Authors show that, compared to the CRS scheme, the RRS scheme can achieve similar performance as long as the resources in the aggregation phase are not very limited.

In its turn, non-orthogonal multiple access (NOMA) has recently attracted a lot of attention as a promising technology for the coming 5G networks to significantly improve the spectral efficiency of mobile communication networks, and/or to meet the demand of massive connectivity. 
Authors in <cit.> present NOMA as a promising solution to support a massive number of MTDs in cellular
networks, while tackling the main practical challenges and future research directions. Actually, NOMA is considered in <cit.> when serving a massive number of devices in cellular-based MTC networks. Energy-efficient clustering and medium access control are investigated in <cit.> to minimize device energy consumption and prolong network battery lifetime, while authors in <cit.> study the throughput and energy efficiency of the NOMA scenario with a random packet arrival model and derive the stability condition for the system to guarantee the performance. 
The key idea of NOMA is to exploit the power domain for multiple access, such that multiple users can be multiplexed at different power levels but at the same time/frequency/code. SIC is utilized to separate superimposed messages at the receiver side <cit.>. 
In general NOMA can be applied on both, downlink and uplink.
A downlink NOMA transmission system is studied in <cit.>, where the authors show that the outage performance of NOMA depends critically on the choices of the targeted data rates and allocated power in scenarios with randomly deployed users; while a dynamic power allocation scheme is proposed in <cit.> to downlink and uplink NOMA scenarios with two users with various quality of service requirements. A NOMA scheme for uplink that allows more than one user to share the same subcarrier without any coding/spreading redundancy is proposed in <cit.>, while the authors establish an upper limit on the number of users per subcarrier to control the receiver complexity. Also in the uplink, authors in <cit.> consider a massive uncoordinated NOMA scheme where devices have strict latency requirements and no retransmission opportunities are available.
Finally, a good overview of NOMA, from its combination with multiple-input multiple-output (MIMO) technologies to cooperative NOMA, as well as the interplay between NOMA and cognitive radio, is provided in <cit.>.


Even though all the recent advances, only few papers focus on evaluating the performance of NOMA by using the stochastic geometry, except for <cit.>. However, the inter-cell interference, which is a pervasive problem in most of the existing wireless networks, is not explicitly considered in <cit.>, and neither do many other works on NOMA.  In contrast, authors in <cit.> do consider the inter-cell interference when evaluating the performance of downlink NOMA on coverage probability and average achievable rate, as well as in <cit.> for downlink and uplink NOMA scenarios. 
In fact, and to the best of our knowledge, <cit.> is the only work analyzing NOMA by means of stochastic geometry in uplink mMTC scenarios for large-scale networks. However, even when authors deal with dense networks, aggregation architectures are not considered, and SIC is assumed perfect in the uplink. This paper aims at filling that gap by characterizing the interference, average success probability and average number of simultaneously served MTDs, for uplink mMTC in a large-scale cellular network system overlaid with data aggregators. We focus on the aggregation phase where the MTDs are allowed to share the same orthogonal channel, while the resource scheduling is implemented at the aggregator side as in <cit.> for OMA setups. The main contributions of this work can be listed as follows:


	
  * We introduce a hybrid access protocol, OMA-NOMA, for the aggregation phase of mMTC systems, while we develop a general analytical framework to investigate its performance in terms of average success probability and average number of simultaneously served MTDs. 	
	
  * We extend the resource scheduling schemes RRS and CRS, proposed in <cit.> to deal with the limited resources, to scenarios where several MTDs are allowed to share the same orthogonal channel. As expected, CRS fits better to our setup since it allows performing adequate power control easily, while improving the system performance. Also, we find the power constraints on the MTDs sharing the same channel in order to attain a fair coexistence of our scheme with purely OMA setups, while power control coefficients are found too, so these MTDs can perform with similar reliability.
	 
	
  * We show that, even when the hybrid scheme could lead to a less reliable system with greater chances of outages per MTD, e.g., due to the additional intra-cluster interference, the number of simultaneous active MTDs could be significantly improved for high access demand scenarios, as long as the success probability does not decrease so much. In that sense, our scheme aims at providing massive connectivity in scenarios with high access demand, which is not covered by traditional OMA setups. Additionally, the CRS scheme requires even a lower average power consumption per orthogonal channel and per MTD, than the OMA setup. 
	
  * We attain approximated, yet accurate, expressions when analyzing the CRS scheme. Compared to the time-consuming Monte-Carlo simulations, even heavier for our hybrid scheme than for the purely OMA setup, our analytical derivations allow for fast computation. Also, imperfection when implementing SIC is incorporated in our proposed analytical framework.




The remainder of the paper is organized as follows. Section <ref> presents the system model and assumptions. Sections <ref> and <ref> discuss the RRS and CRS scheduling schemes for our hybrid access protocol, respectively. Section <ref> studies the power consumption of both schemes for OMA and the hybrid protocol.  Section <ref> presents the numerical results, and finally Section <ref> concludes the paper.

Notation: 𝔼[ · ] denotes expectation, while (A) and  (A|B) are the probability of event A, and (A) conditioned on B, respectively. |S| is the cardinality of set S, ||x|| is the Euclidean norm of vector x and mod(a,b) is the modulo operation. _xF_y is the generalized hyper-geometric function <cit.>, ψ(x) is the digamma function <cit.>, Γ(x) is the gamma function, while Q(a,x) and Beta(x,a,b) are the  regularized incomplete gamma function <cit.> and incomplete beta function <cit.>, respectively. ⌊·⌋ and ⌈·⌉ are the floor and ceiling functions, respectively. 𝐢=√(-1) and Im{z} is the imaginary part of z∈ℂ. f_X(x) and F_X(x) are the Probability Density Function (PDF) and Cumulative Distribution Function (CDF) of random variable (RV) X, respectively. X∼Exp(1) is an exponential distributed RV with unit mean, e.g., f_X(x)=exp(-x) and F_X(x)=1-exp(-x); while Y∼Poiss(m̅) is a Poisson distributed RV with mean m̅, e.g., (Y=y)=1y!m̅^yexp(-m̅) and F_Y(y)=Q(y+1,m̅).



§ SYSTEM MODEL AND ASSUMPTIONS


Consider a cellular network overlaid with data aggregators, which are spatially distributed according to a homogeneous Poisson point process (PPP), denoted Φ_a with density λ_a. Each aggregator, which could function as an ordinary cellular user for certain BS, serves multiple MTDs located nearby. Thus, the result is a cluster point process uniquely defined as:


    Φ⋃_𝐰∈Φ_a𝐰+ℬ^𝐰,


where Φ_a is the parent PPP and ℬ^𝐰 denotes the offspring point process where each point at 𝐬∈ℬ^𝐰 is i.i.d. around the cluster center 𝐰∈Φ_a with distance distribution f_r_a(r_a)=2r_aR_a^2, e.g., uniformly distributed in the disk region of radius R_a. K∼Poiss(m̅) is the instantaneous number of MTDs requiring service in each aggregator, e.g., number of points in ℬ^𝐰, thus, the process is a Matrn cluster point process[A Matrn cluster point process is appropriated when considering either static or low mobility MTDs beign served by the aggregators. Some use cases for mMTC over cellular are: smart utility metering and industry automation <cit.>.]. Notice, that each MTD is associated with a single aggregator even when it could be located within the aggregation areas of several aggregators.


We focus on the uplink where the MTDs across the entire network are served through the same set of orthogonal channels, 𝒩, available at each aggregator, with |𝒩|=N. Differently from <cit.>, here the same orthogonal channel could be used for more than one MTD. Some aggregators could be allocating one MTD per channel because the access demand is not so high, but some other aggregators could be allocating more MTDs per channel to face the increasing access demand, thus we propose and assess a hybrid OMA-NOMA multiple access scenario. Notice that the same aggregator could have channels operating with only one MTD, while others are operating with more. The maximum number of users per orthogonal channel is L, where L=1 reduces to the OMA scenario analyzed in <cit.>, which is used here as a benchmark. We focus on the L=2 setup[This is for analytical tractability, but notice that even when some existing results show that NOMA with more devices may provide a better performance gain <cit.>, this may not be practical. The reason is that considering processing complexity for SIC receivers, especially when SIC error propagation is considered, 2-users NOMA is actually more practical in reality <cit.>.], although some of our results hold for any L≥1. Notice that for L>1 there is both: inter-cluster interference, e.g., interference from MTDs in the serving zone of other aggregators; and intra-cluster interference, e.g., interference from MTDs within the serving area of the same aggregator. After aggregating the MTDs' data, each aggregator relays the entire information to its associated BS. However, the focus of our current  work is on the aggregation phase and we assume that this phase occurs synchronously in all aggregators. Fig.<ref> shows a snapshot of the considered network model. The silent MTDs are those out of the N· L available resources being used by the active MTDs.



We assume the quasi-static fading channel model, where the channel power coefficients, q∈{h,g}, are exponentially distributed with unit mean, e.g., Rayleigh fading. q=h denotes the fading experienced by the signals coming from within the same serving zone, while q=g denotes the fading experienced by the inter-cluster interfering signals. The instantaneous received power at the receiver side is thus given by p_tqr^-α, where p_t is the transmit power from the transmitter, r is the distance between the receiver and transmitter, and α represents the path-loss exponent. Full channel state information (CSI) is assumed at receiver side as in <cit.> in order to obtain benchmark results, and all MTDs are assumed to use statistical full inversion power control <cit.>, which guarantees a uniform user experience while saving valuable energy, with receiver sensitivity ρ. Since here we analyze an interference-limited scenario, ρ does not impact the performance of the network. The aggregators implement the resource scheduling according to one of the schemes presented in the following sections, and the MTDs being considered are those with access granted to the aggregators, since the random access in the network is assumed to be performed[There are some solutions using NOMA for the random access stage as well, e.g., <cit.>. In fact, the work in <cit.> proposes a NOMA scheme where the devices transmit their messages over randomly selected channels, while the random access and data transmissions phases are combined. Readers will realize that our resource scheduling schemes can be easily incorporated to that strategy to improve the overall system performance; however, the details of such implementation are out of the scope of this work.], as in <cit.>.



§ RRS FOR THE HYBRID ACCESS




Herein we explore the RRS scheduling scheme, for which the CSI is only required at the aggregators when decoding the information and not for resource scheduling. Due to its simplicity, this scheme is used as benchmark when compared to a more evolved strategy discussed in Section <ref>, where the benefits of the CSI acquisition are exploited also for resource allocation. In Subsection <ref> we attain the success probability for each MTD in a given channel, which depends on the Laplace transform of the inter-cluster interference. An accurate approximation of the latter is given, which is fundamental to efficiently evaluate the success probability. In Subsection <ref> we characterize the overall system performance.

Under the RRS scheme, N out of the K instantaneous MTDs requiring transmissions are independently and randomly, chosen and matched, one-to-one, with the channels in 𝒩. If K≤N, all MTDs get channel resources, and even N-K channels will be unused. Otherwise, if K>N,  the channel allocation is executed again by allowing the remaining MTDs to share channels with the already served MTDs. This process is executed repeatedly until all the MTDs are allocated or the maximum number of MTDs per channel, L, is reached for all the channels. 


	The Probability Mass Function (PMF) of the number of MTDs allocated to the same channel is given in (<ref>) at the top of the page.
	


	See Appendix <ref>.

The point process of the active MTDs on certain channel n∈𝒩 is obviously a subset of Φ, defined in (<ref>), and can be defined as


    Φ_n⋃_𝐰∈Φ_a𝐰+ℬ^𝐰_n,


where ℬ^𝐰_n∈ℬ^𝐰  denotes the offspring point process with instantaneous number of points around  the cluster center 𝐰∈Φ_a obeying (<ref>). Thus, the generating function of the number of active MTDs on certain channel in one cluster is 


    G_0(z)=∑_u=0^Lc_uz^u,


where c_u=(U=u), while c̅=∑_u=1^Luc_u is the mean number of active MTDs in each channel of certain cluster.



 §.§ MTD success probability for L=2


The intra-cluster interference coming from the MTDs sharing the same channel is faced with SIC[SIC, as in <cit.>, is a common assumption when dealing with NOMA.]. Thus, the SIR[Other than the real SIR at the receiving antenna of an aggregator, we are more interested in the SIR after SIC that can be used to calculate the success probability. Notice that in the L=2 case, the first decoded MTD does not need to perform interference cancellation and directly treats the signal from the second MTD as interference, while the second MTD has to decode first the signal from the other MTD to remove it, which is performed with an efficiency of 100×(1-μ)%.], SIR^r_j,u, of the jth MTD being decoded on a typical channel of certain cluster, given the number of MTDs u sharing the same channel and the RRS scheme is


    SIR^r_1,1   =h/I_r,
    SIR^r_1,2   =max(h',h”)/I_r+min(h',h”),
    SIR^r_2,2   =min(h',h”)/I_r+μmax(h',h”),


where I_r=∑_x∈Φ_n'gr_a^αx^-α is the inter-cluster interference for the RRS scheme with x∈Φ_n'⊂Φ_n denoting both the location and the interfering MTD which occupies certain channel n in other clusters, and r_a is the distance between the MTD and its serving aggregator. Also, μ∈[0,1] is used to model the impact caused by imperfect SIC <cit.>, while h' and h” are the channel power coefficients of both MTDs sharing the channel when u=2. 
Notice that we cannot weight the power of the coexistent nodes on the same channel since CSI information is not exploited for resource scheduling when using the RRS scheme. Also, lim_I_r→ 0SIR_1,2^r is unbounded, but lim_I_r→ 0SIR_2,2^r≤1/μ since min(h',h”)≤max(h',h”), thus the performance of the second MTDs being decoded on the channel is strongly limited by the SIC imperfection parameter.

Assuming a fixed rate coding scheme where the receiver decodes successfully whenever (SIR≥θ), where θ is the SIR threshold, e.g., information rate of log_2(1+θ) (bits/symbol), we state the following theorem.


	The RRS success probability, p^r_j,u, of the jth MTD sharing a typical channel conditioned on u MTDs, is given by
	
	
    p^r_1,1   =ℒ_I_r(θ),
    
    	
    	p^r_1,2   ={[ 2/1+θℒ_I_r(θ)-1-θ/1+θℒ_I_r(2θ/1-θ),                            if 0≤θ<1;                      2/1+θℒ_I_r(θ),                             if θ≥ 1 ].,
    
    	p^r_2,2   ={[ 1-θμ/1+θμℒ_I_r(2θ/1-θμ),                if 0≤θμ<1;                       0,                 if θμ≥ 1 ].,

	where
		
    ℒ_I_r(s)   =exp(2πλ_a∫_0^∞r_𝐰(∑_u=0^Lc_uΥ(r_𝐰,s)^u-1)dr_𝐰),
	
	is the Laplace transform of RV I_r and
	
	
    Υ(r_𝐰,s)   =1/π R_a^2∫_0^R_a∫_0^2πr_adωdr_a/1+sr_a^α(r_𝐰^2+r_a^2+2r_𝐰r_acos(ω))^-α/2.




	See Appendix <ref>.


	Although for some pairs (θ,μ) it is possible to achieve a similar reliability of both MTDs sharing the same channel, e.g., p^r_1,2≈ p^r_2,2, this is not attained in general, and it could be a main drawback when using the RRS scheme and certain homogeneity in the Quality of Service (QoS) is expected. Another interesting fact is that p^r_1,1=p^r_1,2=ℒ_I_r(θ) for θ=1, while p^r_2,2 would be smaller even when μ→ 0, e.g., ℒ_I_r(2θ).
	Additionally, notice that μ<1/θ is required such that the second user being decoded on certain channel has chances to succeed. Thus, the greater the SIR threshold, θ, the greater the impact of imperfect SIC, μ.



	Expression (<ref>) for L≥ 2 is upper (almost surely) and lower bounded by 
	
    ℒ_I_r(s)   a.s≤ℒ_I_r^up(s)=exp(-χ s^2/α∑_u=1^Lc_uu^2/α),
    ℒ_I_r(s)   ≥ℒ_I_r^lo(s)=exp(-χc̅ s^2/α),

	where χ=1/2λ_aπ R_a^2Γ(1+2α)Γ(1-2α), while 
	
    ℒ_I_r(s)   ≈β_0ℒ_I_r^up(s)+β_1ℒ_I_r^lo(s),

	
	provides an approximation with β_0,β_1∈[0,1] and β_0+β_1=1. 

 
 
 	See Appendix <ref> for derivation of (<ref>) and (<ref>), while (<ref>) is a weighted average of both bounds, thus attaining an approximation that will be more accurate than at least one of the bounds. The weighted average becomes relevant if we know a priori which of the bounds in (<ref>) and (<ref>) is more accurate for the setup being analyzed and we give it a greater weight, e.g., β_0≶β_1, otherwise trivial choice β_0=β_1=0.5 is advisable, as we adopt here.
 
 
 
 
 Fig. <ref> shows the bounds and approximation attained in Theorem <ref> while comparing them with the exact value given in (<ref>).  Specifically, Fig. <ref> shows the performance for α∈{3,5}, while we set c_2=1, which conduces to scenarios where the bounds in (<ref>) and (<ref>) are the least tight since T_1 in (<ref>) has no impact on (<ref>), e.g., T_1=1. As noticed, the greater the α the lesser the tightness of the bounds, however the approximation performs very well whatever the setup, which is also appreciated in Fig. <ref>.  Notice that for relatively small N, the approximation in (<ref>) is accurate even when λ_a increases and the bounds are not tight. When N increases, c_0 and c_1 increase, which increases the contribution of T_1 in (<ref>), thus increasing the tightness of the bounds, and the accuracy of the approximation. 
 


 §.§ Overall Performance for L=2


The following two lemmas characterize the system performance in terms of average over all MTD success probability, and average number of simultaneously served MTDs since one of the main benefits of NOMA techniques is the possibility of offering service to a great number of devices simultaneously.


	Conditioned on being using certain channel, the average over all MTD success probability is
	
	
    p^r_succ=c_1/1-c_0p^r_1,1+c_2/2(1-c_0)(p^r_1,2+p^r_2,2).

	



	Knowing that c_u=(U=u) it is straightforward to attain (<ref>).





	The average number of simultaneously served MTDs is given in (<ref>) at the top of the page, where
	
	
    A_1   =m̅Q(N+1,m̅)-exp(-m̅)m̅^N+1/N!,
        
        A_2   =exp(-m̅)(m̅^2N/(2N-1)!-m̅^N+1/N!)-(m̅-2N)(Q(2N,m̅)-Q(N+1,m̅)),
        
    	A_3   =N(Q(N+1,m̅)+1)

	



	In (<ref>), (a) comes from averaging the number of simultaneously active and successful MTDs, while (b) comes from regrouping terms and using the expressions of the CDF and expected K on one interval (see Appendix <ref>) along with some algebraic transformations.


	When m̅ grows above 2N, the last term in equality (a) of (<ref>) becomes the largest contributor for K̅_r, with K̅_r→ N(p^r_1,2+p^r_2,2)=2Np^r_succ since ∑_k=2N^∞(K=k)≈ 1 and c_2≈ 1 (c_0≈ c_1≈ 0) in (<ref>). Thus, when comparing with an L=1 setup under the same circumstances, e.g., with average number of simultaneously served MTDs being Np_succ^OMA, the condition required for the hybrid access scheme[In general, the required condition can be extended for any L, thus p_succ^HYB>1Lp_succ^OMA.] to overcome the OMA system (L=1) is that p^r_succ>12p_succ^OMA. 
    Notice that the hybrid scheme with RRS could lead to a less reliable system with greater chances of outages per MTD, e.g., due to a larger inter-cluster interference  and additional intra-cluster interference. However, the number of simultaneous active MTDs could be improved for high access demand scenarios, e.g., m̅∼ 2N, as long as the success probability does not decrease so much.	




§ CRS FOR THE HYBRID ACCESS


Herein we explore the CRS scheduling scheme, which contrary to the RRS scheme, strongly relies on the CSI for resource scheduling. Thus, the CRS scheme is more adjusted to NOMA scenarios, where CSI is a keystone for efficiently decoding multiple user data over the same orthogonal channel with SIC <cit.>. Similar than in the previous section, we attain first in Subsection <ref> the success probability for each MTD in a given channel along with necessary approximations. Later, in Subsection <ref>, we find the power constraints on the MTDs sharing the same channel in order to attain a fair coexistence of our scheme with purely OMA setups, while power control coefficients are found too, so these MTDs can perform with similar reliability. Finally, in Subsection <ref> we characterize the overall system performance.

Under the CRS scheme, the MTD with better fading (equivalently, better SIR) will be preferentially assigned with the available channel resources. An aggregator with K instantaneous  MTDs requiring transmission has the knowledge of their fading gains. Let {h_1,...,h_i,...,h_K} denote the decreasing ordered channel gains, where h_i-1>h_i. If K≤ N all the MTDs will be chosen, but if K>N the aggregator will pick the N MTDs with better channel gains, i.e., h_1,...,h_N, and then will assign the channel set 𝒩 to them <cit.>. As a continuation, the remaining MTDs can  be still allocated sharing those same resources, i.e., users N+1,...,K go to the second round for allocation. This process is executed repeatedly until all the MTDs are allocated or the maximum number of MTDs per channel, L, is reached. 
Both Lemma <ref> and the point process of the active MTDs on certain channel n∈𝒩 characterized through (<ref>) and (<ref>) hold here.



 §.§ MTD success probability for L=2


Under the CRS scheme and using SIC to face the intra-cluster interference, the SIR, SIR^c (i)_j,u, of the jth MTD being decoded on a typical channel, given the first MTD allocated there has the ith larger channel coefficient, h_i, and there are u MTDs sharing that same channel, is given by


    SIR^c (i)_1,1   =h_i/I_c,
    SIR^c (i)_1,2   =a_ih_i/I_c+b_ih_i+N,
    SIR^c (i)_2,2   =b_ih_i+N/I_c+μ a_ih_i.


Notice that differently from the RRS scheme, here we can weight the power of coexistent nodes on the same channel through a_i and b_i since CSI is used for resource allocation. Of course, some kind of feedback from the aggregators would be required.
Once again, the performance of the first MTD being dedoced in the channel is somewhat unbounded, e.g., lim_I_c→ 0SIR_1,2^c (i) is unbounded, while the performance of the second one is not,  but this time we can relax that situation by choosing a_i<b_i since lim_I_c→ 0SIR_2,2^c (i)<b_i/μ a_i.
However, weighting the transmit power of the MTDs sharing the same channel conduces to a marked point process with marks m_x=1 when the MTD on x is alone in the channel, which happens with probability c_1, and m_x∈{a_i,b_i} for those MTDs sharing the same channel, which happens with probability c_2. Thus, the inter-cluster interference for the CRS scheme is I_c=∑_x∈Φ_n'gm_xr_a^αx^-α.
By letting a_i+b_i=δ be a fixed value we impose some kind of total transmission power constraint <cit.>. This is crucial for NOMA scenarios, and here it is particular important in order to control the interference from the inter/intra-cluster MTDs sharing the same channel.  
Now, assuming that the receiver can decode successfully (SIR exceeds a threshold θ), we state the following theorem.


	The CRS success probability, p^c (i)_j,u, of the jth MTD being decoded on a typical channel, given the first MTD allocated there has the ith larger channel coefficient, h_i, and there are u MTDs sharing that same channel, is approximately given by
	
	
    p^c (i)_j,u   ≈1/2-1/π∫_0^∞1/φIm{ℒ_I_c(-𝐢φ)exp(-𝐢φ B_j,u^(i,K))}dφ,

	
	where
	
	
    ℒ_I_c(s)   =exp(2πλ_a∫_0^δ∫_0^∞r_𝐰(c_0+c_1Υ(r_𝐰,s)+c_2Υ(r_𝐰,a_is)Υ(r_𝐰,(δ-a_i)s)-1)f_a_i(a_i)dr_𝐰da_i),
       ≈∑_t∈{1,2}β_t-1exp(-χ(c_1+c_2t^α-2/αδ^2/α)s^2/α),

	is the Laplace transform of RV I_c for L=2, Υ(r_𝐰,s) is defined in (<ref>), and β_0,β_1∈[0,1], β_0+β_1=1, such that (<ref>) with β_0=1 and β_1=1 provide, almost surely, upper and lower bounds for (<ref>), respectively. Finally,

	
    B_1,1^(i,K)   =ψ(K+1)-ψ(i)/θ,
    
    	B_1,2^(i,K)   =(a_i/θ-b_i)ψ(K+1)+b_iψ(i+N)-a_i/θψ(i),
    
    	B_2,2^(i,K)   =(b_i/θ-μ a_i)ψ(K+1)+μ a_iψ(i)-b_i/θψ(i+N).




	See Appendix <ref>.

Notice that even in the case when all MTDs operating on the same channel in the network are using the same a_i and b_i coefficients, e.g., a_i and b_i with deterministic values, evaluating (<ref>) using (<ref>) is not an easy task. A closer look at that expression makes us suspect on its efficiency. This is because it requires evaluating two inner triple integrals.
In fact, several numerical tests we run corroborated that evaluating (<ref>), using (<ref>) with fixed a_i,b_i values, is highly inefficient and computationally too heavy. Thus, the approximation given in (<ref>)  becomes necessary, which holds also under the premise that the pair (a_i,b_i) has not to be simultaneously the same for all MTDs operating on the same channel. This would be required for instance in order to optimize in some way the system performance. Additionally, notice that B_1,2^(i,K) and B_2,2^(i,K) defined in (<ref>) and (<ref>), respectively, are functions of the power control coefficients a_i and b_i of the typical links, which are assumed fixed for the entire period of evaluation.




Fig. <ref> shows the accuracy of (<ref>) by plotting the exact values using (<ref>) for fixed a_i,b_i and δ=1 and comparing also with simulations by choosing a_i, b_i uniformly random. The values for a_i, b_i are interchangeable since the perceived interference is the same. The remaining values for the system parameters are the same than the previously used when discussing Fig. <ref>. Both, Fig. <ref> and Fig. <ref> showing ℒ_I_c as a function of s and N, respectively, corroborate the idea behind  (<ref>). This is, the system performance depends heavily on δ rather than the individuals a_i, b_i, or their distribution. Notice that the exact values for two completely different pairs (a_i,b_i) are very similar, and even when a_i, b_i were chosen randomly, the performance is kept alike. 



	We know that the interference under the CRS scheme with a_i=b_i=1 is characterized in the same way as in the RRS case, ℒ_I_c=ℒ_I_r, since no weights to the transmit power, e.g., no marks, are assigned; and notice that (<ref>) captures well this phenomena since with δ=2 we attain (<ref>). This means that whatever the values of a_i, b_i, as long as δ=2 the interference is distributed approximately equal as for the RRS scenario.



	An alternative, and easy of evaluating, expression for the CRS success probability, p_j,u^c (i), for a given δ, is
	
	
    p_j,u^c (i)≈1/2-∑_t∈{1,2}β_t-1/π∫_0^∞exp(-ς_tφ^2/α)sin(ϱ_tφ^2/α-φ B_j,u^i,K)/φdφ,
	
	
	where ς_r=ν_rcos(πα), ϱ_r=ν_rsin(πα) and ν_r=χ(c_1+c_2r^α-2/αδ^2/α).
	


	The idea here is substituting into (<ref>) the approximation for ℒ_I_c(s) given in (<ref>). Using (-𝐢)^2/α=cos(πα)-𝐢sin(πα) and Im{pexp(-q𝐢)}=-psin(q), along with simple algebraic transformations, we attain (<ref>). The summation comes from both sum terms in (<ref>).


	


	The CRS success probability, p_j,u^c, of the jth MTD sharing a typical channel conditioned on u MTDs, is given in (<ref>) and (<ref>) at the top of the next page, where j∈{1,2}.



	Note that (a) comes directly from 
	the total probability theorem, and also averaging p_j,u^c (i) over all possible values of i given k and N. On the other hand, (b) comes from using the expressions for (K≤ N) and (K=k) easily obtained through the PDF and CDF of K. Also, (k,N)=k-N if N+1 ≤ k≤ 2N-1, and the approximation is because we substituted the infinite sum  for a finite sum until k_max to reach (<ref>).


Notice that a proper value for k_max is such that the sum over the remaining k>k_max does not contribute significantly to the success value in (<ref>). Thus, we could choose k_max such ∑_k=k_max+1^∞(K=k)<τ→ Q(k_max+1,m̅)>1-τ holds, where τ is the maximum allowable error we admit when approximating with (<ref>), e.g., τ=10^-5 and m̅=30→ k_max=56.



 §.§ Practical issues


	If we look closely at (<ref>) and (<ref>) we can notice there must exist some choice for a_i and b_i in order to attain a similar reliability for both MTDs sharing the orthogonal channel. Thus, we resort to SIR_1,2^c (i)=SIR_2,2^c (i) with a_i+b_i=δ. However, the knowledge of the inter-cluster interference is required and it is a major drawback for this 

		method. Alternatively, we state the following theorem avoiding that problem.

	


	A proper approximate choice for a_i and b_i in order to attain a similar reliability for both MTDs sharing the channel when u=2 is given by
	
	
    a_i   =δ(1+1/θ)(ψ(K+1)-ψ(i+N))/(1+μ+2/θ)ψ(K+1)-(μ+1/θ)ψ(i)-(1+1/θ)ψ(i+N),
    
    	b_i   =δ-a_i.



	Notice that p_1,2^c (i) and p_2,2^c (i) in (<ref>) only differ in the terms B_1,2^(i,K) and B_2,2^(i,K). Thus, solving B_1,2^(i,K)=B_2,2^(i,K) with a_i+b_i=δ conduces to an approximate choice for these parameters. 



	The values of a_i and b_i aside of the index i, strongly rely on the instantaneous number of MTDs contesting for transmission resources, K, the number of available channels, N, the SIR threshold, θ, and the SIC imperfection parameter, μ. These values for a_i, b_i do not guarantee the same instantaneous SIR  for both MTDs sharing the channel, but when averaging over a long period[Assuming that each of them will be occupying the same ith channel and the same order when transmitting in future rounds.] it guarantees a similar reliability for them. Also, all MTDs have the same chances of occupying any of the channels in 𝒩, thus for high loaded systems where m̅>N, e.g., the probabilities of all the channels being occupied by two MTDs are great, an average similar reliability will be attained. 




On the other hand, both the RRS scheme, and CRS scheme with relatively large δ, do not favor the coexistence with purely OMA setups. The reason is because these NOMA setups would increase the interference, e.g., up to twice and approximately δ larger for RRS and CRS schemes, respectively, caused to the OMA setups, and since this is not compensated by multiplexing several users, the performance of the OMA clusters is affected. On the other hand, even when only our hybrid OMA-NOMA scheme is utilized, it is expected that all aggregators will not be under the same average access demand, e.g., same m̅.
	Therefore, those with low demand could be operating with OMA almost all the time and the larger interference of previously schemes will be impractical. 
Thus, limiting the interference is crucial and that could be done by properly selecting a relatively small δ.



	The required δ, δ^*, for a fair coexistence[We refer to “fair coexistence” when for any aggregator, the interference coming from the outside topology (the inter-cluster interference) remains the same regardless of the alternative (OMA protocol or our hybrid approach) the outside clusters utilize.] between OMA and our hybrid setup is approximated by the solution of	
	
    ξ^δ^2/α-1+ξ^2^α-2/αδ^2/α-1=2,

	where ξ=exp(-χ c_2 s^2/α), and it is bounded by
	
	
    2^2-α/2≤δ^*≤ 1.




	See Appendix <ref>. 


	Interestingly, even though δ^* depends completely on the system parameters e.g., λ_a, R_a, c_u, θ and α, we were able to limit its range only as a function of α. The greater the α the smaller δ^* and consequently the greater the limitation on NOMA setup over nodes. Also, the fact that δ^*< 1 is expected, means that the NOMA setup has to operate with an overall consumption power inferior to OMA setup. 




 §.§ Overall Performance for L=2


Differently than in the RRS case, where we were able to attain the average over all MTD success probability based only on a linear combination of p_j,u^r as stated in (<ref>), here we cannot do the same since the weights of p_1,1^c and p_j,2^c when averaging are different for each k and dependent on the index of the ordered channels, i. Instead, the average over all MTD success probability is given in (<ref>) as a function of p_j,u^c (i)  at the top of the page, where (a) comes from using the expressions for (K≤ N), (K=0) and (K=k).

On the other hand, the average number of simultaneous served MTDs is given in (<ref>), shown below of (<ref>) on top of the page, where A_1 was given in (<ref>) and the approximation is because the same reasons than previously discussed. Notice that Remark <ref> holds here for the CRS scheme as well.





§ POWER CONSUMPTION ANALYSIS



	The average transmit power per orthogonal channel for the OMA (L=1) and, RRS and CRS of our hybrid scheme (L=2) is given by
	
	
    p̅_t   ={[                      (1-c_0)Ψ,                        for OMA;                   (c_1+δ c_2)Ψ for our hybrid scheme with L=2 ].,

	where Ψ=2ρ R_a^α/α+2.



	For OMA, RRS and CRS with L=2, the average transmit power is 𝔼[(1-c_0)ρ r^α]=(1-c_0)ρ𝔼[r^α], 𝔼[c_1ρ r^α+2c_2ρ r^α]=c̅ρ𝔼[r^α] and 𝔼[c_1ρ r^α+c_2(a_i+b_i)ρ r^α]=(c_1+δ c_2)ρ𝔼[r^α], respectively. For RRS c̅=c_1+δ c_2 since δ=2, and now it is only necessary to compute 𝔼[r^α],
	
	
    𝔼[r^α]=∫_0^R_ar^αf_r(r)dr=2/R_a^2∫_0^R_ar^α+1dr=2R_a^α/α+2=Ψ/ρ,

	
	and (<ref>) is attained.


Obviously, the RRS scheme with L=2 will always require a greater power consumption than an OMA setup, since whenever two MTDs are transmitting on the same channel of a cluster, the power consumption doubles. On the other hand, the CRS scheme allows to reduce the power consumption by adopting a relatively small δ, e.g., δ^*. This behavior is illustrated in Fig. <ref> for c_0=0. All the schemes have a common starting point in Ψ since all of them are equivalent when no MTDs require to share the same orthogonal channel, e.g., c_2=0→ c_1=1. When c_2 increases, the power consumption of the RRS and CRS schemes more, as long as δ 2, reaching the maximum difference for c_2=1. Notice that choosing δ^* we guarantee operating with the same interference than an OMA setup while the power consumption reduces since δ^*<1.






§ NUMERICAL RESULTS


Both, simulation and analytical results, are presented in this section in order to investigate the performance of our hybrid scheme as a function of the system parameters while comparing it with an OMA setup. The analytical results for the RRS scheme come from using the exact expressions; while for the CRS scheme we use the approximations. Unless stated otherwise, results are obtained by setting m̅=60, λ_a=10^-4.4/m^2 (39.81/km^2), R_a=40m, α=3.6, μ=0 and θ=1. The value of δ^* is found by solving numerically (<ref>) whenever is required. Simulation results are generated using 50000 Monte Carlo runs and a sufficiently large area such 400 aggregators are placed on average.


Fig. <ref> shows the success probability of the MTDs sharing the same channel conditioned on u=2, e.g., (<ref>) and (<ref>) for RRS, and (<ref>) for CRS. The idea is to show how fair the schemes are when allocating the transmission resources. Notice that for the RRS scheme, the gap between both MTDs performance keeps quite constant, independently of N. This is because that gap relies strongly on the gap between max(h',h”) and min(h',h”), (see (<ref>) and (<ref>))[Notice the gap also relies strongly on the value of μ.], which only depends on the number of MTDs requiring transmissions and not on the available channels. While for the CRS schemes with fixed a_i, b_i, the gap tends to widen when increasing N because h_i becomes larger with respect to h_i+N in (<ref>) and (<ref>). Of course, a given system is projected to work given one value of N, and by properly choosing some fixed a_i and b_i we can reduce the performance gap between the MTDs sharing the channel, which is not possible for the RRS scheme since weighting the transmit powers is not available. Notice also that by choosing a_i and b_i according to (<ref>) and (<ref>) both MTDs reach similar performance, hence the fairest scheme. 



Fig. <ref>a shows the average over all MTD success probability for the different schemes and L={1,2}. Each scheme with L=1 performs better than for L=2. Also, as expected and in general, the CRS setups outperforms the RRS scheme, except when choosing a_i, b_i according to (<ref>) and (<ref>), for which a greater fairness is attained instead. Notice that when N is small, CRS outperforms RRS; while as N increases, their curves tend to overlap. This is due to the fact that when N is large,

	i) e.g, greater than m̅ for L=1, most of the time the number of MTDs is less than the available resources such that the implementation of the CRS is almost the same as the RRS; 

ii) e.g., greater than m̅/2 for L=2, the impact of ordering the channels becomes less significant since most of the time all the MTDs will be allocated.  While when N increases more and more, the system tends to behave as when L=1.

The change in the CRS curves from decreasing to increasing  occurring close to m̅/2 is somewhat explained by that latter argument. Even when achieving a higher reliability with an L=2 setup is not possible, we are able to enhance the average number of simultaneously served MTDs, as shown in Fig. <ref>b. The success probability does not deteriorate and a significant improvement on K̅ is attained, fundamentally when N is not large. Notice that this advantage is evinced when L=2, especially with the CRS scheme, so to cover a high instantaneous access demand m̅>N, which is even more favorable than predicted in Remark <ref>. By setting L=2, the number of efficiently served MTDs could be up to twice the number with a single MTD per orthogonal channel setup. Particularly attractive are the configurations operating with δ^* since no change in the perceived interference occurs when switching transmission schemes from L=1 and L=2 setups and vice versa.



Fig. <ref> shows the average number of simultaneously served MTDs as a function of the density of aggregators for μ∈{0,10%}, Fig. <ref>a, and the SIC imperfection coefficient, Fig. <ref>b. When both, the network is sparse and the SIC imperfection coefficient are not so restrictive, the performance of the L=2 setup increases. Since N=30, each channel per cluster is operating with two MTDs almost all the time, which are more sensitive to the interference, hence their performance will be affected if either λ_a, θ, or even R_a are larger. If c_2 was smaller, e.g., larger N, this situation becomes less critical, although the gap between the L=1 and L=2 setups would be smaller as shown previously in Fig. <ref>. Notice that the imperfect SIC degrades the system performance when L=2, but even for a high imperfection such 10%, the advantage over the OMA setup keeps evident for a wide range of values for the system parameters. For instance, the CRS scheme with a_i=b_i=δ^*/2 overcomes the CRS scheme with L=1 when λ_a≲ 1.3· 10^-4 for μ=0, while λ_a≲ 1· 10^-4 would be required for μ=10%. Since SIC is only related with the L=2 setup, the OMA setup appears shown as a straight line for both RRS and CRS schemes in Fig. <ref>b. Notice that the setup with fixed power coefficients, e.g., RRS and CRS with fixed a_i, b_i, are the most affected when μ increases since those coefficients work well for certain system parameters but others will be required if they change, e.g. different μ in this case. It is expected that a smaller a_i, hence larger b_i, work better as μ increases (see (<ref>) and (<ref>)). The setup with a_i and b_i given by (<ref>) and (<ref>), respectively, adapts better to the different μs since this is a parameter taken into account when performing their calculations. This is, the values of a_i and b_i that allow a similar reliability between both MTDs sharing the same channel in a given cluster, depend on μ, thus no additional adjustment on them is required. It is clear that failing to efficiently  eliminate the intra-cluster interference could reduce significantly the benefits from NOMA, and can be a  challenging issue for implementing NOMA in practice.


Finally, notice that simulations and analytical expressions, even when several of them are approximations, fit well in all the cases, e.g., Fig. <ref>-<ref>, which validates our findings[Analytical expressions, even when some of them are numerically evaluated, have an additional value since simulations of these scenarios require a huge amount of computational resources specially for L>1.].



§ CONCLUSION


In this paper, we analyzed the uplink mMTC in a large-scale cellular network system overlaid with data aggregators. We propose a hybrid access scheme, OMA-NOMA, while developing a general analytical framework to investigate its performance in terms of average success probability and average number of simultaneously served MTDs for two scheduling schemes RRS and CRS. Also, we found the power constraints on the MTDs sharing the same channel to attain a fair coexistence with purely OMA setups, while power control coefficients are found too, so that both MTDs can perform with similar reliability.  Our analytical derivations allow for fast computation compared to the time-consuming Monte-Carlo simulations, which are even heavier for our hybrid scheme than for a purely OMA setup. The numerical results show that


	
  * our hybrid access scheme aims at providing massive connectivity in scenarios with high access demand, which is not covered by traditional OMA setups, and even with lower average power consumption per orthogonal channel and per MTD, the hybrid scheme with CRS outperforms the OMA setup; 
	
  * failing to efficiently  eliminate the intra-cluster interference could reduce significantly the benefits from NOMA, and can be a  challenging issue for implementing NOMA in practice;
	
  * our mathematical derivations, besides being easy to evaluate, are accurate.	


Future work could focus on the relaying phase, while finding some strategies to cope with the larger aggregated data. Additionally, we intend to deeply investigate strategies in order to optimally decide when to switch from a purely OMA setup to our hybrid scheme.




§ PROOF OF LEMMA <REF>

The PMF of the number of MTDs sharing the same channel, u, conditioned on the number of MTDs requiring transmissions, k, is given by


    (U=u|k)=   {[                          1,           for u=L if k≥ N L;                 1-k/N+⌊kN⌋,         for u=⌊kN⌋ if k<N L;                   k/N-⌊kN⌋, for  u=⌊kN⌋+1=⌈kN⌉ if k<N L;                          0,                   otherwise ].
    .


Notice that if k is equal or greater than the number of available resources, NL, there will be for sure u=L MTDS per channel in the representative cluster. Otherwise, if k<NL, only two consecutive values for u are possible. For instance, if N=10, L=4 and k=26, then 2 MTDs will be allocated in each of 4 channels (2× 4=8 MTDs), and 3 MTDs in the remaining 6 channels (3× 6=18 MTDs). Therefore, the probability of one channel being occupied by 2 MTDs is 4/10=0.4=1-26/10+⌊ 26/10 ⌋, while the probability of being occupied by 3 is the complement, 0.6=26/10-⌊ 26/10 ⌋.

Now, the required PMF, (U=u)=(U=u|k)(K=k), can be calculated as follows

    (U=u)   ={[               ∑_k=NL^∞(K=k),                      for u=L; (1-k/N+⌊kN⌋)∑_k=0^NL-1(K=k),                   for u=⌊kN⌋;   (k/N-⌊kN⌋)∑_k=0^NL-1(K=k),           for  u=⌊kN⌋+1=⌈kN⌉;                           0,                    otherwise ],
    .
       ={[                             ∑_k=0^N-1(1-kN)(K=k),                    u=0; ∑_k=Nu^N(u+1)-1(1-kN+u)(K=k)+∑_k=N(u-1)^Nu-1(kN-u+1)(K=k),   u=1,...,L-1;                          ∑_k=NL^∞(K=k)+∑_k=N(L-1)^NL-1(kN-L+1)(K=k), u=L;               0,                                               otherwise ].
    .


Using the PDF of K, (K=k), its CDF, (K≤ k), and the expected K on one interval,

    ∑_k=a^k=bk(K=k)=m̅(Q(a,m̅)+Q(1+b,m̅))+exp(-m̅)(m̅^a/(a-1)!-m̅^1+b/b!)
    
       a=0=m̅Q(1+b,m̅)-exp(-m̅)m̅^1+b/b!b→∞=m̅(1-Q(a,m̅))+exp(-m̅)m̅^a/(a-1)!,

and regrouping similar terms, we attain (<ref>). 


§ PROOF OF THEOREM <REF>


We proceed as follows


    p^r_1,1   =(SIR_1,1>θ)=𝔼_I[(h>θ I|I)]=𝔼_I[exp(-θ I)|I],
    
    p^r_1,2   =(SIR_1,2>θ)
    =𝔼_I[(max(h',h”)-θmin(h',h”)>θ I|I)]
    =𝔼_I[(v_1>θ I|I)]
       (a)={[ 𝔼[2/1+θexp(-θ I)-1-θ/1+θexp(-2θ/1-θ I)|I],                                   if 0≤θ<1;                       𝔼[2/1+θexp(-θ I)|I],                                    if θ≥ 1 ].,



    p^r_2,2   =(SIR_2,2>θ)
    =𝔼_I[(min(h',h”)-θμmax(h',h”)>θ I|I)]
    =𝔼_I[(v_2>θ I|I)]
       (b)={[ 𝔼[1-θμ/1+θμexp(-2θ/1-θμI)|I],                     if 0≤θμ<1;                            0,                      if θμ≥ 1 ].,


where v_1=max(h',h”)-θmin(h',h”) and v_2=min(h',h”)-θμmax(h',h”), while (a) and (b) come from using their CDF expressions, which are given by


    F_V_1(v_1)   ={[ 1-
    2/1+θexp(-v_1)+1-θ/1+θexp(-2v_1/1-θ),                                 if 0≤θ<1;                        1-2/1+θexp(-v_1),                                  if θ≥ 1;                                          ].,
    
    F_V_2(v_2)   ={[ 1-
    1-θμ/1+θμexp(-2v_2/1-θμ),                    if 0≤θμ<1;                           1,                     if θμ≥ 1;                              ].,

for v_1,v_2>0. Notice that the last equalities in (<ref>), (<ref>) and (<ref>) are equivalent to (<ref>), (<ref>) and (<ref>), respectively.

For the derivation of (<ref>) we have to use the fact that the Poisson cluster process defined in (<ref>) is a Neyman-Scott process <cit.> with Probability Generating Functional (PGFL) <cit.>
	
	
    G[υ]=exp(λ_a∫_ℝ^2(G_0(∫_ℝ^2v(𝐰+y)f_y(y)dy)-1)d𝐰).
	

Based on g∼Exp(1) which allows to state


    v(𝐰+y)   =𝔼_g[exp(-sgr_a^α||𝐰+y||^-α)]=1/1+sr_a^α||𝐰+y||^-α,

and according to the cosine law, e.g., ||𝐰+y||=(r_𝐰^2+r_a^2-2r_𝐰r_acos(ω))^1/2, while substituting (<ref>) and PDFs expressions f_r_a(r_a)=2r_aR_a^2 and f_ω(ω)=12π into (<ref>), we reach ℒ_I_r(s) in (<ref>). 



§ PROOF OF THEOREM <REF>


Performing some algebraic transformations on (<ref>) we have


    ℒ_I_r(s)   =exp(2πλ_a∫_0^∞r_𝐰(∑_u=0^Lc_uΥ(r_𝐰,s)^u-1)dr_𝐰)
    
       (a)=exp(2πλ_a(∫_0^∞r_𝐰(∑_u=0^L(c_uΥ(r_𝐰,s)^u-c_u))dr_𝐰))
       (b)=exp(2πλ_a∑_u=1^Lc_u(∫_0^∞r_𝐰(Υ(r_𝐰,s)^u-1)dr_𝐰))
       (c)=exp(2π c_1λ_a∫_0^∞r_𝐰(Υ(r_𝐰,s)-1)dr_𝐰)_T_1exp(2πλ_a∑_u=2^Lc_u(∫_0^∞r_𝐰(Υ(r_𝐰,s)^u-1)dr_𝐰))_T_2,


where (a) comes from ∑_u=0^Lc_u=1, (b) is attained by regrouping terms, and (c) by pulling out the term associated with u=1. Notice that T_1 matches the Laplace transform of an HPPP with density c_1λ_a and one active MTD per channel per cluster, which is given by <cit.>


    T_1=exp(-χ c_1s^2/α)


On the other hand, T_2 includes the contribution of the clustered MTDs, e.g., u≥ 2. For each u, the related term in T_2 matches the Laplace transform of an HPPP with density c_uλ_a and u active MTDs per channel per cluster. It has been observed in <cit.> that the SIR complementary CDFs, e.g., Laplace transform of the interference, for different point processes appear to be merely horizontally shifted versions of each other (in dB), as long as their diversity gain is the same. Thus, scaling the threshold s by this SIR gain factor (or shift in dB) G, we have


    ℒ_I(s)≈ℒ_I,ref(s/G).


However, G is also a function of s but for many setups it keeps approximately constant and consequently it can be determined by finding its value for an arbitrary value of s <cit.>.
Using the PPP as the reference model, the limit of G as s→ 0, G_0, is relatively easy to calculate <cit.>


    G_0=MISR_PPP/MISR,
 

where the MISR is the mean of the interference-to-(average)-signal ratio IS̅R=I𝔼_h[S]. Since S=h here and 𝔼[h]=1, we have that MISR=𝔼[I]. Now, considering the contribution of each u in T_2 separately we have


    G_0,u=𝔼[I_PPP]/𝔼[I_r]=1/u.
 

Notice that 𝔼[I_PPP] and 𝔼[I_r] are divergent measures for our system model because the resultant point process would no be locally finite since we assumed a path loss model ||x||^-α <cit.>. However, the quotient depends merely on the density of both process and since λ_PPP=c_uλ_a and λ=uc_uλ_a we attain the last equality in (<ref>). Obviously, ℒ_I,PPP(us) works almost surely[Except for relatively large s (<ref>) might serve as an upper bound but it becomes an accurate approximation for (<ref>).] as an upper bound for ℒ_I_r(s) in (<ref>) if G_0<1 (see <cit.> for a geometrical perspective), which holds here since u≥ 2.

  Now, using <cit.> as the HPPP of reference with one fixed MTD per orthogonal channel[The clustered process with one MTD per orthogonal channel (per cluster) is indeed a HPPP with density λ_a because the displacement theorem in stochastic geometry<cit.>.], we attain 
  
  
    T_2a.s≤exp(-χ∑_u=2^Lc_u (us)^2/α)

  
  with asymptotic equality as s→ 0. Substituting (<ref>) and (<ref>) into (<ref>) yields (<ref>). 

On the other hand, the lower bound comes from the corresponding HPPP with same intensity that for each u in T_2, e.g., uc_uλ_a. This is because in our system model the performance in terms of success probability of an HPPP with the same intensity that any clustered process will be always worse since we are expecting a greater number of close interfering nodes. Once again using <cit.>, we have


    T_2≥exp(-χ∑_u=2^Luc_us^2/α).


Now, substituting (<ref>) and (<ref>) into (<ref>) yields (<ref>). 

For both (<ref>) and (<ref>) the equality fits to the asymptotic cases λ_a→∞,0, R_a→∞,0.  


§ PROOF OF THEOREM <REF>


Let's assume the case where K>N since for K≤ N the system behaves exactly as in the RRS scheme and the problem is already solved, e.g., p^c (i)_j,u=p^r_j,u. According to high order statistic theory <cit.>, the PDF of the ith best channel power gain is


    f_H_i,K(x)   =K!exp(-ix)(1-exp(-x))^K-i/(i-1)!(K-i)!.


Now we have


p^c (i)_j,u=(SNR^c (i)_j,u>θ)=(Θ^i_j,u>θ I_c),


where Θ^i_1,1=h_i, Θ^i_1,2=a_ih_i-θ b_ih_i+N and Θ^i_2,2=b_ih_i+N-θμ a_ih_i. 

Unfortunately, the distribution of Θ^i_j,u, even for the simplest case Θ^i_1,1, conduces to very complicated expression for the CDF, preventing us to follow the same path we used for the RRS scheme. Instead we proceed as follows
	
	
    p^c (i)_j,u   =(I_c<Θ^i_j,u/θ)(a)=1/2-1/π∫_0^∞1/φIm{ℒ_I_c(-𝐢φ)𝔼_Θ^i_j,u[exp(-𝐢φΘ^i_j,u/θ)]}dφ
       (b)≈1/2-1/π∫_0^∞1/φIm{ℒ_I_c(-𝐢φ)exp(-𝐢φ𝔼[Θ^i_j,u]/θ)}dφ,

	
	where (a) comes from the Gil-Pelaez inversion theorem <cit.> and the approximation in (b) comes from the Jensen inequality. Also,
	
	
    𝔼[Θ^i_1,1]   =𝔼[h_i]=∫_0^∞xf_H_i,K(x)dx
    	=K!/(i-1)!(K-i)!∫_0^∞xexp(-ix)(1-exp(-x))^K-idx
       =(-1)^K-iK!/K^2(i-1)!(K-i)!(K^2xBeta(exp(x),-K,K-i+1)
                    -exp(-Kx) _3F_2(-K,-K,i-K,1-K,1-K,exp(x)))|_x=0^x=∞
       =Γ(i)Γ(K-i+1)(ψ(K+1)-ψ(i))/(i-1)!(K-i)!=ψ(K+1)-ψ(i),

	
	while it is straightforward obtaining 𝔼[Θ^i_1,2] and 𝔼[Θ^i_2,2] from (<ref>), yielding (<ref>)-(<ref>). Substituting them into (<ref>) we attain (<ref>). 
	

	To attain (<ref>) we require to use the same arguments than previously discussed when deriving the result in (<ref>). However, now the point process has marks 1,a_i,b_i, with b_i=δ-a_i, and we require to include the marks along with their probabilities when evaluating (<ref>) <cit.>.   
	
		Notice that it is intractable evaluating (<ref>) efficiently when a_i,b_i are random since requires an additional integration[For instance, the uniform distribution, which is very simple and probably unrealistic for the scenario discussed, with PDF f_a_i(a_i)=1δ, is already cumbersome when evaluating (<ref>).] in the exponent, hence we propose using some approximations as follows		

    ℒ_I_c(s)   (a)≈exp(2πλ_a∫_0^∞r_𝐰(c_0+c_1Υ(r_𝐰,s)+c_2/2(Υ(r_𝐰,a_is)^2+Υ(r_𝐰,b_is)^2)-1)dr_𝐰)
       (b)=exp(2π c_1λ_a∫_0^∞(Υ(r_𝐰,s)-1)dr_𝐰)exp(2πc_2/2λ_a∫_0^∞(Υ(r_𝐰,a_is)^2-1)dr_𝐰)exp(2πc_2/2λ_a∫_0^∞(Υ(r_𝐰,b_is)^2-1)dr_𝐰),


where (a) comes from using fixed values of a_i,b_i, thus, avoiding the additional integration, and using the relation between the geometric and arithmetic mean, x+y2≥√(xy)→ xy≤x^2+y^22, as an approximation, with equality when a_i=b_i; while (b) comes from regrouping terms and ∑_u=0^2c_u=1. Now, by using the same procedure than previously discussed in the proof of Theorem <ref>, e.g., finding upper and lower bounds for (<ref>), and averaging both, 
we attain


    ℒ_I_c(s)≈∑_t∈{1,2}β_t-1exp(-χ(c_1+c_2(t/2)^α-2/α(a_i^2/α+b_i^2/α))s^2/α),


where β_0,β_1∈[0,1] and β_0+β_1=1.
Since 2α<1, and the relation between the generalized mean (or power mean) with exponents 2α and 1, the following result holds


    (a_i^2/α+b_i^2/α/2)^α/2≤a_i+b_i/2
    
    a_i^2/α+b_i^2/α≤2(δ/2)^2/α.

By substituting (<ref>) into (<ref>) we attain (<ref>). See Fig. <ref> for more insights on the accuracy. 



§ PROOF OF THEOREM <REF>
	
In order to attain a fair coexistence between the OMA and NOMA setups the interference must be kept the same. Thus, we need to match (<ref>) with the Laplace transform of the interference for an equivalent OMA setup as shown next.


    ∑_t∈{1,2}β_t-1exp(-χ(c_1+c_2t^α-2/αδ^2/α)s^2/α)   =exp(-χ(1-c_0)s^2/α)
    exp(-χ c_1s^2/α)(exp(-χ c_2 (sδ)^2/α)+exp(-χ c_2 2^α-2/2(sδ)^2/α))   (a)=2exp(-χ(1-c_0)s^2/α)
    exp(-χ c_2 (sδ)^2/α)+exp(-χ c_2 2^α-2/2(sδ)^2/α)   (b)=2exp(-χ(1-c_0-c_1)s^2/α)
    ξ^δ^2/α+ξ^2^α-2/αδ^2/α   (c)=2ξ


where (a) comes from setting β_0=β_1=0.5 and evaluating the left-hand sum,  (b) comes from dividing both terms by exp(-χ c_1 s^2/α), and (c) by 1-c_0-c_1=c_2 and setting ξ=exp(-χ c_2 s^2/α). Dividing both terms by ξ we attain (<ref>). The solution is an approximation since (<ref>) is an approximation. 
Now, bounding the solution is simple since (<ref>) is a combination of upper, t=1, and lower, t=2, bounds. Thus,



			exp(-χ						(c_1+c_2δ^2/α)s^2/α)   ≥exp(-χ(1-c_0)s^2/α)

			c_1+c_2δ^2/α   ≤1-c_0

			c_2δ^2/α   ≤ c_2
δ   ≤ 1,    
			exp(-χ (c_1+c_2 2^α-2/αδ^2/α)s^2/α)   ≤exp(-χ(1-c_0)s^2/α)

    		c_1+c_2 2^α-2/αδ^2/α   ≥ 1-c_0

    		2^α-2/αδ^2/α   ≥ 1
δ   ≥ 2^2-α/2,


completing the proof. 


IEEEtran


