I Introduction
Wireless powered backscatter communications (WPBC) have been widely adopted in ultralowpower communication, such as radiofrequency identification (RFID) tags and Internet of Things (IoT) [duan2020ambient]. In general, a WPBC system consists of three components: a carrier emitter, a tag reader, and multiple backscatter tags. At first, the carrier emitter generates radio frequency (RF) signals to activate backscatter tags. Then, each tag harvests the energy of incident RF signals for its internal circuit operations and modulates its information bits by simply reflecting the received carrier signals with different antenna impedances. Finally, the tag reader receives the backscattered signals and decodes them. Thanks to the low power consumption of tags, WPBC enables operation in a batteryfree manner. Moreover, instead of generating RF signals by a tag itself, a backscatter tag does not require conventional RF components like local oscillators, mixers, and converters, thereby significantly facilitating the implementation in IoT.
Ia Related Works and Motivation
Despite the idea of backscatter communications was conceived by Stockman as early as in 1948 [Stockman1948], it is not put into practice until recent years. In 2012, a practical backscatter receiver equipped with a nonlinear nearoptimal detector was designed in [kimionis2012bistatic], and subsequently the same authors proposed to deploy multiple carrier emitters to increase the coverage of wireless sensor networks [kimionis2014increased]. Furthermore, a backscatter network was designed in [liu2017full] to enable both oneway wireless power transfer (WPT) and twoway wireless information transfer. These works are devoted to WPBC applications in either pointtopoint scenario or smallscale networks.
The communication range of WPBC is very limited due to the small amount of energy harvested by a backscatter tag [MXia2015TSP]. To extend the coverage of WPBC, dedicated power beacons can be employed. By using a stochastic geometry approach, the works [han2017wirelessly] and [bacha2018backscatter] studied the coverage and capacity of largescale WPBC networks. In [zhu2018inference]
, an intelligent backscatter sensor system was developed, where machine learning techniques were leveraged to enhance signal processing of backscatter sensors. In
[lu2018wirelesspowered], a hybrid devicetodevice (D2D) system integrating ambient backscattering with wirelessly powered communications was introduced, where two mode selection protocols were devised so as to adapt to diverse propagation environments. However, although power beacons benefit shortening the distance of WPT, they increase substantially the infrastructure investment of network operators [Tharindu2018Simultaneous].Thanks to its low cost and mobility, unmanned ground vehicle (UGV) is an attractive medium to replace power beacons in WPBC for largescale IoT networks. In [Shuai2017Wirelessly], a UGV was used to act as a mobile relay in wireless powered twoway communication system with two terminals. It is demonstrated that, with an appropriate trajectory design, a UGV can enlarge the achievable data rate region of the system. In [Shuai2019Backscatter], the optimal trajectory planning and power allocation at a UGV in a backscatter system with multiple terminals was investigated, to achieve an energy balance between UGV motion and data transmission. To the best of our knowledge, the trajectory planning of UGV and corresponding resource allocation for largescale networks with massive terminals is still an open problem. To resolve this problem, in this paper, we design a UGVassisted backscatter system for largescale IoT networks, where after careful network planning, a spiral trajectory for the UGV is specified, as well as several efficient resourceallocation methods.
IB Main Contributions
This paper develops a novel network architecture where a UGV is employed to enable a largescale backscatter network. Specifically, an AP serving as both data transmission/reception and energy carrier emitter is deployed at the center of the network while a backscatter reader is mounted on a UGV that can move along a predesigned trajectory. With the motion of UGV, the reader can visit all tags in sequence, then collect the backscattered data, and finally send them back to the AP for further processing. To fully exploit massive antennas at the AP, both halfduplex (HD) and fullduplex (FD) operation modes are studied and compared. In summary, the main contributions of this paper are as follows:

A UGVassisted WBPC architecture suitable for largescale IoT networks is devised. In the network, the AP transmits RF signals to activate tags. Then, the UGV moves along a predesigned trajectory and the reader mounted on the UGV collects the backscattered data. Finally, the reader relays data to the AP for further processing;

To optimize the trajectory of UGV, the network coverage is tessellated into hexagonal cells and the optimal radius of cells is derived analytically under both FD and HD modes. Also, the optimal network planning under FD mode is proved identical to that under HD mode;

Under HD mode, the optimal resource allocation of the network is performed, including Tx/Rx beamforming at the AP and dynamic power allocation at the tag reader;

Under FD mode, a successive convex approximation (SCA) based algorithm is designed to obtain the optimal power allocation at the tag reader and Tx/Rx beamforming at the AP. Moreover, two suboptimal schemes with low complexity or in closed form are developed.

Simulation results show that dynamic power allocation at the tag reader is much more energyefficient than Tx beamforming at the AP. Also, the AP operating in FD mode has higher energy efficiency than that in HD mode.
The proposed system together with the developed algorithms can be applied, e.g., in vehicle manufacturing, where massive tags are deployed to monitor the status of various industrial equipments and UGVs are already onsite.
To detail the contributions described above, the remainder of this paper is organized as follows. Section II develops the network architecture, including network modeling and signal modeling. Section III formulates the problem of the joint design of network planning and resource allocation, given that the AP operates in HD mode. Then, Section IV investigates the case that the AP operates in FD mode, where selfinterference is accounted for. Section V presents and discusses simulation results and, finally, Section VI concludes the paper.
Notation
: Vectors and matrices are denoted by lower and uppercase letters in boldface, respectively. Calligraphic letters indicate sets or optimization problems, depending on the context. The operators
and represent the floor and ceiling functions of a real number , respectively, while returns the real part of a complex number . The operators , and take the absolute value, Euclidean norm and expectation of , respectively. The superscripts and mean the inverse and Hermitian transpose of a matrix, respectively. The symbols and indicate the real and complex spaces with dimension , respectively. The abbreviation implies thatfollows a circularly symmetric complex Gaussian distribution with mean
and variance
. Finally, andindicate an identity matrix with proper size and the set of positive integers, respectively.
Ii The Proposed Network Architecture
This section starts with network modeling, followed by the corresponding signal modeling.
Iia Network Modeling
We consider a large number of backscatter tags uniformly distributed in coverage area
with density . An AP serving as a central processing unit or an edge computing node is located at the center of the coverage area. Each tag has no builtin battery and can only communicate with a nearby tag reader by backscattering the signals transmitted by the AP. To serve such a largescale network, a tag reader is mounted on a UGV and shares its battery with the UGV, and the UGV moves along a carefully designed trajectory and visits all tags in succession. During each sojourn time of the UGV, the reader communicates with its nearby tags and then relays the collected data to the AP. In other words, the tag reader serves as a mobile decodeandforward relay.As shown in Fig. 1, by using a similar idea to cellular communications, we tessellate the network coverage area into hexagonal cells. Suppose that the AP is located in the central cell with ID , then, around the central cell, there are regular hexagonal cells and they are indexed in a spiral form with cell ID , as shown by the red curve with arrows in Fig. 1. For ease of notation, is also used to denote the center of the cell as well. Moreover, all the cells are grouped into layers, indexed with from the inner to outer layer. For instance, the six grey cells around the central cell form the first layer (i.e., ) while the twelve yellow cells adjacent to the first layer constitute the second layer (i.e., ). Clearly, for the layer, there are cells in total.
Due to the relatively short distance between those tags in the central cell and the AP, they are allowed to communicate directly with the AP, or, for security purposes, the central cell can be seen as an exclusion zone without any tag [MXia2015TSP]. As a result, the UGV starts to move from the cell with ID , and then to the cell with ID , until the last cell with ID in sequence (cf. Fig. 1). At each cell, the UGV sojourns for a while and the reader communicates with the tags in the cell and, then, relays the collected data to the AP.
It is assumed that the reader and tags are all equipped with a single antenna while the AP has antennas. The propagation channels between the reader and the AP are assumed subject to Rayleigh fading and the largescale pathloss has a pathloss exponent . In this section, we study the case that the AP operates in halfduplex mode, namely, all antennas at the AP are used for either data transmission or reception in timedivision multiplexing. For comparison purposes, the fullduplex mode will be investigated in Section IV.
Figure 2 illustrates the frame structure and workflow of the network shown in Fig. 1, where the strategy of timedivision multiple access (TDMA) is adopted among the tags. Specifically, as shown in Fig. 2(a), a time block of length is divided into slots and the slot corresponds to the sojourn duration of the UGV in the cell, for all . Each slot is further divided into subslots and the subslot is devoted to the communication between the AP and the tag in the cell, where . In particular, during the subslot when the UGV sojourns in the cell, that is, in the subslot, the communications between the AP, the tag, and the reader consists of three phases:

Downlink Phase: The AP serving as carrier emitter broadcasts signal and the tag harvests energy from its received signal;

Backscattering Phase: The tag modulates its data on the same RF carrier as the AP and backscatters them to the reader;

Uplink Phase: The reader decodes the backscattered signal and relays the decoded data to the AP for further processing.
For a given network and UGV trajectory described above, it is clear that the energy consumed by the motion of UGV depends on the number of cells or, equivalently, the radius of hexagonal cells. For an energyconstrained IoT network, the radius of hexagonal cells becomes the dominant factor of the network performance. Before embarking on the technical details of determining cell radius, insofar we assume that the cell radius is a prior, say, , then, the area of each hexagonal cell can be given by . Accordingly, the relationship between the coverage area , the number of cells , the number of cell layers , and the radius can be explicitly established as
(1) 
Assuming that the tags are uniformly distributed in the coverage area with density , the average number of tags in each cell is
(2) 
By accounting for the motion time and the sojourn time of the UGV spent on all cells, the total number of subslots that the UGV completes visiting all tags in the network, i.e., the length of time block , is explicitly computed as
(3) 
where denotes the motion speed of the UGV, in the unit of meter per subslot.
Now, we compute the energy consumption of the network. On the one hand, in light of (1)(3), the energy consumed by the UGV can be expressed as
(4) 
where and are parameters used in the mobility model of a UGV, e.g., for a Pioneer 3DX robot [mei2006deployment]; denotes the Tx power of tag reader for relaying data back to the AP when it sojourns at the cell, and refers to the power consumed by internal circuits of both the UGV and tag reader [Shuai2016Multipair]. Clearly, the first term on the righthand side (RHS) of (4) denotes the motion energy of UGV; the second term indicates the Tx energy of tag reader, and the last term accounts for the energy consumed by internal circuits of both the UGV and tag reader. On the other hand, let denote the beamforming vector of AP intended for the tag in the cell, then, the total energy consumed by the AP can be expressed as
(5) 
where the first term on the RHS of (5) indicates the Tx energy of AP, and in the second term stands for the power consumed by internal circuits of AP.
As illustrated in Fig. 2(b), the horizontal distance between the AP and the center of the cell can be easily calculated as , where is the layer number in which the cell is located. Given the height of the AP, the TxRx distance between the AP and UGV is given by
(6) 
Furthermore, thanks to the relatively small radius of the cells, the distance between the AP and the tag in the cell, is approximated by the distance between the AP and the center of the cell, that is, , . To make the approximation as accurate as possible, a tolerable pathloss difference in dB is specified as
(7) 
It is noteworthy that the tolerance plays a pivotal role in the following network planning since the values of and depend on the cell radius .
IiB Signal Modeling
As described above, in each subslot the communications between the AP, the tag, and the reader consists of three transmission phases, namely, the downlink phase, the backscattering phase, and the uplink phase. Now, we elaborate the corresponding signal model at each phase.
IiB1 Downlink Phase
Let with normalized energy (i.e., ) be the symbol transmitted by the AP to the tag. Then, the received signal at the tag can be written as
(8) 
where denotes the complex channel vector between the AP and tag , and indicates the circular sysmmetric additive white Gaussian noise (AWGN) at tag .
IiB2 Backscattering Phase
When the UGV sojourns at cell , the reader mounted on the UGV receives the backscattered signals from tag as well as the interfering signals from the AP. Let with denote the Tx symbol of tag , then, under ideal synchronization, the received signal at the reader can be expressed as
(9)  
(10) 
where denotes the power reflection coefficient of the tag; is the complex channel coefficient between the tag and tag reader; denotes the complex channel vector between the AP and tag reader, and is an AWGN at the reader. Notice that, the approximation in (10) is introduced by neglecting the backscattered noise , which is plausible due to its relatively lower power compared to the local noise [Fuschini2008Analytical]. Moreover, due to the relatively short distance between the reader and tags compared with that between the reader and AP, the pathloss between the reader and tags is not accounted for.
Since the backscattered signals suffer from both downlink attenuation and backscatter attenuation, at the tag reader the strength of interfering signals from the AP is generally much stronger than that of the backscattered signals. Specifically, in light of (10), the stronger signal, i.e., , can be decoded firstly and then subtracted from the Rx signal [bharadia2015backfi, gong2018backscatter]. Thus, the resultant signal can be written as
(11) 
Accordingly, the Rx signaltonoiseratio (SNR) at the reader is given by
(12) 
IiB3 Uplink Phase
Let the Tx signal of the reader be with , which contains the decoded information from tag in the previous phase. Then, the received signal at the AP from tag is
(13) 
where and denote the Tx power of the tag reader and the Rx beamforming vector of the AP, respectively; is an AWGN at the AP. As a result, the received SNR at the AP can be computed as
(14) 
In light of (12) and (14), the achievable data rate between the AP and tag can be explicitly given by
(15) 
where the factor before the logarithm operator is due to the dualhop transmission, namely, from the tag to the tag reader and then to the AP.
Remark 1 (On the CSI acquisition, control signaling, and synchronization).
In the network under study, while the tags are batteryfree and cannot generate active radio signals, the channel state information (CSI) can be acquired at the tag reader through semiblind channel estimation
[Zhang2019A] or blind channel estimation [8320359]. As for the control signaling, following the similar procedure developed in [Moeinfar2012Design], the reader allocates an identification (ID) number to each tag, and then, if necessary, wakes up the desired one by its unique ID and communicates with it by using the RF signals from the AP. As passive tags are usually equipped with lowpower oscillators, in practice, the EPC Gen2 protocol [EPC2018] uses slotted Aloha to synchronize the AP and all tags in a cell. Namely, each tag adjusts its clock offset upon reception from the AP and introduces a guard time to ensure that it communicates with the reader at the allocated subslot [Yildirim2019On]. In practice, even imperfect synchronization has little effect on the network performance [Wang2012Efficient]. For more details on these engineering issues, the interested reader is referred to, e.g., [Wang2012Efficient, 8737551] and references therein.Iii Joint Network Planning and Resource Allocation
In this section, we investigate the network planning and optimal resource allocation at the AP and tag reader, including Tx/Rx beamforming at the AP and Tx power allocation of the tag reader. As the data transmission of tag reader consumes much less energy than the motion of UGV, the objective of our design is to minimize the total energy consumption of the AP and UGV while satisfying the qualityofservice (QoS) requirement of tags and the AP. Accordingly, let , , and , the optimization problem of joint network planning and resource allocation can be formulated as follows:
(16a)  
(16b)  
(16c)  
(16d)  
(16e)  
(16f)  
(16g)  
(16h) 
where and in (16a) are earlier defined in (4) and (5), respectively; in (16b), in (16c), and in (16d) indicate the minimum required data rate for tag , the maximum amount of energy stored in the UGV, and the maximum allowable system operating time, respectively; in (16f) and in (16g) represent the maximum Tx power of the reader and of the AP, respectively; finally, (16h) indicates that Rx beamforming vectors are normalized in power.
Since (16e) puts constraint on the cell radius (i.e., ), which further influences the number of cells (i.e., ), the number of constraints of depends heavily on the tolerable pathloss difference . As a result, is hard to solve, if not impossible. To tackle this difficulty, we first express in terms of using (1):
(17) 
then, inserting (17) into (3) yields
(18) 
where denotes the motion time of UGV, which is a function of the number of cell layers (i.e., ); and where stands for the sojourn time of UGV, which depends only on the number of tags in the network (by recalling the fact that timedivision multiplexing is applied to tags) but is independent of ; and the approximation is introduced to reach (18), which holds if . Consequently, substituting (18) into (4) and (5) as well as performing some algebraic manipulations, Problem can be equivalently rewritten as
(19a)  
(19b)  
(19c)  
(19d)  
While Problem is still nontrivial as most variables are coupled together, after careful observation, we find that the number of cell layers (i.e., ) can be decoupled from others and determined separately, as detailed below.
Iiia Optimal Number of Cell Layers
As illustrated in Fig. 1, more cells introduce longer trajectory of the UGV and, thus, higher motion energy consumption. Since data transmission consumes much less energy than the motion of UGV, the energy consumed by both the UGV and tag reader is dominated by the motion of UGV, or equivalently by the number of cells. From this perspective, the number of cells should be as few as possible so as to save energy (strictly speaking, is an increasing function of , as proved in Appendix A). On the other hand, as shown in (16e), the cell radius is limited by the tolerable pathloss difference. As a result, the optimal radius of hexagonal cells is the maximal value that satisfies the equality of (16e) and the optimal number of cell layers can be determined, as summarized in the following theorem.
Theorem 1.
Proof.
See Appendix A. ∎
With the resultant , the optimal radius of hexagonal cells (i.e., ), the number of cells (i.e., ), and the average number of tags in each cell (i.e., ) can be explicitly determined as
(21)  
(22)  
(23) 
Once the optimal number of layers (i.e., ) is determined according to Theorem 1, the cell radius (i.e., ) can be easily computed through (21) and, finally, the spiral trajectory of the UGV is fixed as per Fig. 1. This completes the task of optimal network planning. Next, we investigate the dynamic resource allocation.
To avoid excessive notation, in the rest of this section we remain to use and , instead of and with the understanding that their optimal values are derived. After is determined by (20), Problem can be simplified as
(24a)  
(24b)  
where denotes the amount of energy used for data transmission. It is clear that except (16b), Problem is decoupled in and , both in the objective function and the constraints. To break the coupling of (16b), on account of (IIB3), it is straightforward that (16b) is equivalent to
(25)  
(26) 
where and . Therefore, can be decomposed into two independent subproblems, one for and the other for :
(27)  
(28) 
In the following, is first solved to obtain the Tx beamforming vector and, then, is solved to obtain both the Rx beamforming vector and the power allocation vector .
IiiB Optimal Tx Beamforming of the AP
Since all tags work in a timedivision multiplexing fashion, can be further decomposed into subproblems, one for each :
(29a)  
(29b)  
(29c) 
Thanks to the special structure of , can be analytically derived, as formalized below.
Theorem 2.
The optimal Tx beamforming vector of the AP for the tag is given by
(30) 
Otherwise, Problem is infeasible.
Proof.
See Appendix B. ∎
It is observed that (30) is essentially the wellknown conjugate beamforming. Also, the condition in (30) implies that, if the channel quality between the AP and the tag is too poor to satisfy the QoS requirement under the peak Tx power constraint of AP, as well as
is infeasible, and the tag is not allowed to make transmission at the moment.
IiiC Optimal Rx Beamforming of the AP and Optimal Tx Power of the Tag Reader
As the objective function of shown in (28) is monotonically increasing with , the term in (26) must be maximized so as to minimize the Tx power of tag reader (i.e., ) while satisfying (26). Together with (16h), the design of is the maximal ratio combining (MRC) problem [roy2004maximal] while satisfying . Consequently, the optimal Rx beamforming vector of AP in slot is given by
(31) 
With determined by (31), can be simplified as
(32a)  
(32b)  
Since smaller would make (16f) and (24b) easier to satisfy and lead to small objective function of Problem , it is straightforward that the optimal should satisfy the equality in (32b), as summarized in the following theorem:
Theorem 3.
The optimal Tx power of the tag reader in cell is given by
(33) 
if
(34) 
Otherwise, Problem is infeasiable.
The first condition in (34) implies that, if the channel quality between the tag reader and AP is too poor to satisfy the QoS requirement under the reader’s peak Tx power constraint, the reader shall suspend its communication to the AP at the moment. The second condition in (34) means that the total energy used for data transmission in the network must be no larger than the prepared energy , defined immediately after (24). If one or both of the conditions cannot be satisfied, as well as is infeasible.
Iv FullDuplex Mode
Insofar we have assumed an AP with antennas operating in HD mode. Compared to HD mode, FD multiplexing achieves higher spectral efficiency and lower latency [Gang2014In]. In particular, as illustrated in Fig. 3(a), the antennas at the AP are divided into two parts: antennas are used for data transmission and the remaining ones are for data reception, with . As the Tx signal of the AP is generally much stronger than its Rx signal, the FD mode suffers from severe selfinterference [wang2020performance]. While some advanced signal processing techniques can be exploited to suppress selfinterference, there is still residual selfinterference. To account for this residual selfinterference in our design, a selfinterference channel between the Tx and Rx antennas is modeled as , which is invariant to different time slots since the antennas are relatively fixed [everett2014passive, shende2017half]. In the following, we investigate the optimal network planning and resource allocation where the AP operates in FD mode while suffering from residual selfinterference. Unlike the HD mode, to deal with the residual selfinterference inherent in FD mode, the Rx beamforming vectors at the AP must be carefully redesigned.
Due to the selfinterference, the signal model in FD mode differs significantly from the preceding HD counterpart. More specifically, in the uplink phase the Rx signal at the AP for the tag is now reexpressed as
(35) 
where the second term in the parentheses denotes the selfinterference caused by the FD operation of the AP. Then, the instantaneous Rx signaltointerferenceplusnoiseratio (SINR) at the AP can be computed as
(36) 
and accordingly, the achievable data rate between the AP and the tag is modified from (IIB3) as
(37) 
Therefore, the optimization problem of joint network planning and resource allocation under FD mode is the same as that in , except that in (16b) and in (16h) are replaced with and , respectively.
As previously stated in Subsection IIIA, since the data transmission of tag reader consumes much less energy than the motion of UGV, even though the channel conditions under FD mode differ from those under HD mode, the network planning can follow the same procedure as what we have done in the preceding section. In other words, under FD mode, the radius of hexagonal cells, the number of cells, and the number of tags in each cell are also determined by (21), (22), and (23), respectively. Moreover, for notational simplicity, we remain to use and , instead of and , to denote the number of cells and average number of tags in each cell, respectively. Then, similar to the derivation from to , the new problem introduced in this section can also be derived as with and being replaced, which is now labeled as .
In light of given by (IV), the discontinuous constraint (16b) is equivalent to the following two inequalities:
(38)  
(39) 
where and . Although the discontinuity is resolved, the problem is still challenging to solve since the variables , , and are coupled all together in (39).
To proceed, we observe that is only involved in (39), regardless of either the objective function or other constraints of . Accordingly, to minimize the transmit power in (39), can be determined as per
(40) 
It is not hard to recognize that (40) is in the form of generalized Rayleigh quotient [horn2012matrix, Sec. 4.2]. Thus, the solution to (40) can be explicitly expressed as
(41) 
Inserting (41) into (39) and after performing some straightforward manipulations, (39) can be rewritten as
(42) 
After eliminating the nontrivial matrix inversion using ShermanMorrison formula [hager1989updating], (42) can be equivalently rewritten as:
(43) 
where the term is convex (for more details, please refer to Appendix C). Finally, the optimization problem can be reduced w.r.t. only and :
(44)  
The problem is still nonconvex due to the concave constraint (38) and the concave term in (IV). In principle, by using the technique of semidefinite relaxation (SDR) [Mohammadi2016Throughput, song2018joint, chalise2019full], the nonconvex can be transformed into a quasiconvex optimization problem, and then be solved by the bisection method [boyd2004convex]. However, as the SDRbased method requires lifting the problem to higher dimensional space, it is computationally inefficient. Even worse, the SDRbased method may fail when the AP’s antenna size becomes larger. To tackle the ineffectiveness of SDRbased method, we use the SCAbased method to reformulate , as detailed below.
Iva SCAbased Joint Optimization (JOSCA)
For the convexification of (38), the term is expanded as the firstorder Taylor series. Then, (38) at the iteration can be rewritten in terms of the optimal value of at the previous iteration (i.e., ) as
(45) 
which is linear. Likewise, (IV) can be rewritten as in the following convex form:
(46) 
where is the firstorder Taylor series of the term defined in (IV). Therefore, to solve reduces to solve a sequence of convex problems, where at the iteration the problem can be written as
(47)  
which can be easily solved by using popular convex optimization package like CVX [boyd2004convex]. With the optimal solutions to , for all , we set and then is solved to produce . This process iterates until it converges. In summary, based on the SCA, the overall procedure for solving the original problem is formalized in Algorithm 1.
Computational Complexity: Taking CVX for instance to solve , at each iteration, using the interiorpoint method yields complexity [ben2001lectures, Sec. 6.6.3]. Thus, the overall complexity of Algorithm 1 is , where denotes the number of iterations.
IvB Suboptimal Method with Equal Power Allocation (SOEPA) at the Reader
While the SCAbased Algorithm 1 is more efficient than the SDRbased method, it becomes computationally challenging when the size of the problem continues to grow. For instance, considering an IoT network for vehicle manufacturing, there must be a huge number of tags or sensors for equipment and environment monitoring, asset location, and inventory control, among others. As explicitly computed above, the computational complexity of Algorithm 1 is proportional to with being the average number of tags in each cell. This high complexity becomes impractical if is large. To suit massiveaccess networks, in the following, we develop a suboptimal parallel algorithm whose complexity is independent of .