Backscatter Data Collection with Unmanned Ground Vehicle: Mobility Management and Power Allocation
Shuai Wang Minghua Xia Yik-Chung Wu



[4]















theoremTheorem
theorembox 200 
lemmaLemma
lemmabox 500 
propositionProposition
propositionbox 500 
corollaryCorollary
corollarybox 500 
propertyProperty
propertybox 500 
remarkRemark
remarkbox 500 
claimClaim
claimbox 500 







Backscatter Data Collection with Unmanned Ground Vehicle: Mobility Management
and Power Allocation
    Shuai Wang, Minghua Xia, and Yik-Chung Wu      

This work was supported by the National Natural Science Foundation of China (NSFC) under Grant No. 61671488, and by the Major Science and Technology Special Project for “New Generation Communication and Network”, Guangdong Province, China.

S. Wang and Y.-C. Wu are with the Department of Electrical and Electronic Engineering, The University of Hong Kong, Hong Kong (e-mail: {swang, ycwu}@eee.hku.hk).

M. Xia is with the School of Electronics and Information Technology, Sun Yat-sen University, Guangzhou, 510006, China (e-mail: xiamingh@mail.sysu.edu.cn).


    December 30, 2023
=====================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




Collecting data from massive Internet of Things (IoT) devices is a challenging task, since communication circuits are power-demanding while energy supply at IoT devices is limited.
To overcome this challenge, backscatter communication emerges as a promising solution as it eliminates radio frequency components in IoT devices. Unfortunately, the transmission range of backscatter communication is short.
To facilitate backscatter communication, this work proposes to integrate unmanned ground vehicle (UGV) with backscatter data collection.
With such a scheme, the UGV could improve the communication quality by approaching various IoT devices.
However, moving also costs energy consumption and a fundamental question is: what is the right balance between spending energy on moving versus on communication?
To answer this question, this paper studies energy minimization under a joint graph mobility and backscatter communication model.
With the joint model, the mobility management and power allocation problem unfortunately involves nonlinear coupling between discrete variables brought by mobility and continuous variables brought by communication.
Despite the optimization challenges, an algorithm that theoretically achieves the minimum energy consumption is derived, and it leads to automatic trade-off between spending energy on moving versus on communication in the UGV backscatter system.
Simulation results show that if the noise power is small (e.g., ≤-100 dBm), the UGV should collect the data with small movements. However, if the noise power is increased to a larger value (e.g., -60 dBm), the UGV should spend more motion energy to get closer to IoT users.



Backscatter communication, Internet of Things (IoT), mixed integer optimization, quality-of-service (QoS), unmanned ground vehicle (UGV).




§ INTRODUCTION


With a wide range of commercial and industrial applications, Internet of Things (IoT) market is continuously growing <cit.>, and the number of inter-connected IoT devices is expected to exceed 20 billion by 2020.
However, these massive IoT devices (e.g., sensors and tags) are usually limited in size and energy supply <cit.>, making data collection challenging in IoT systems.
To this end, backscatter communication is a promising solution, because it eliminates radio frequency (RF) components in IoT devices <cit.>.
Unfortunately, due to the round-trip path-loss, the transmission range of backscatter communication is limited <cit.>.
This can be seen from a recent prototype in <cit.>, where the wirelessly powered backscatter communication only supports a range of 1 meter at the data-rate of 1 kbps.


To combat the short communication range, this paper investigates a viable solution that the backscatter RF transmitter and tag reader are mounted on an unmanned ground vehicle (UGV).
With such a scheme, the UGV could vary its location for wireless data collection, thus having the flexibility of being close to different IoT devices at different times <cit.>.
However, since moving the UGV would consume motion energy, an improperly chosen path might lead to excessive movement, thus offseting the benefit brought by movement <cit.>.
Therefore, the key is to balance the trade-off between spending energy on moving versus on communication, which unfortunately cannot be handled by traditional vehicle routing algorithms <cit.>, since they do not take the communication power and quality-of-service (QoS) into account.


In view of the apparent research gap, this paper proposes an algorithm that leads to automatic trade-off in spending energy on moving versus on communication.
In particular, the proposed algorithm is obtained by integrating the graph mobility model and the backscatter communication model.
With the proposed model, the joint mobility management and power allocation problem is formulated as a QoS constrained energy minimization problem.
Nonetheless, such a problem turns out to be a mixed integer nonlinear programming problem (MINLP), which is nontrivial to solve due to the nonlinear coupling between discrete variables brought by moving and continuous variables brought by communication.
This is in contrast to unmanned aerial vehicle (UAV) based systems in which only continuous variables are involved <cit.>.
To this end, the optimality condition of the MINLP is first established, which helps in reducing the problem dimension.
Then, an efficient algorithm, which is guaranteed to obtain the global optimal solution, is proposed.
By adopting the proposed algorithm, minimum energy consumption is achieved at the UGV, and simulation results are presented to further demonstrate the performance of the proposed algorithm.

The rest of this paper is organized as follows. In Section II, the system model, which includes the mobility model and the backscatter communication model, is described.
Then, the joint mobility management and power allocation problem is formulated in Section III.
The algorithm for computing the optimal solution is derived in Section IV, and an efficient initialization is proposed in Section V.
Finally, numerical results are presented in Section VI, and conclusions are drawn in Section VII.


Notation.
Italic letters, simple bold letters, and capital bold letters represent scalars, vectors, and matrices, respectively.
Curlicue letters represent sets and |·| is the cardinality of a set.
We use (a_1,a_2,⋯) to represent a sequence and [a_1,a_2,⋯]^T to represent a column vector, with (·)^T being the transpose operator.
The operators Tr(·) and (·)^-1 take the trace and the inverse of a matrix, respectively.
Finally, 𝔼(·) represents the expectation of a random variable.




§ SYSTEM MODEL




 §.§ Mobility Model




We consider a wireless data collection system, which consists of K IoT users and one UGV equipped with a RF transmitter and a tag reader.
The environment in which the UGV operates in is described by a directed graph (𝒱,ℰ) as shown in Fig. 1, where 𝒱={1,⋯,M} is the set of M vertices representing the possible stopping points, and ℰ is the set of directed edges representing the allowed movement paths <cit.>.
To quantify the path length, a matrix 𝐃=[D_1,1,⋯,D_1,M;⋯;D_M,1,⋯,D_M,M]∈ℝ^M× M_+ is defined, with the element D_m,j representing the distance from vertex m to vertex j (D_m,m=0 for any m).
If there is no allowed path from vertex m to vertex j, we set D_m,j=+∞ <cit.>.
Notice that 𝐃 is a constant matrix since the locations of all vertices are pre-determined.


To model the movement of the UGV, we define a visiting path 𝒬=(y_1,y_2,⋯,y_Q) where y_j ∈𝒱 for j=1,⋯,Q and (y_j,y_j+1)∈ℰ for j=1,⋯,Q-1, with Q-1 being the number of steps to be taken.
Without loss of generality, we assume the following two conditions hold:


  (i) y_1=y_Q. This is generally true as a typical UGV management scenario is to have the UGV standing by at the starting point (e.g., for charging and maintenance services) after the data collection task  <cit.>.
    For notational simplicity, it is assumed that vertex y_1=y_Q=1 is the start and end point of the path to be designed.


  (ii) There are no repeating vertices among (y_1,⋯,y_Q-1). This is true because if a vertex m is visited twice, we can always introduce an auxiliary vertex with D_M+1,j=D_m,j and D_j,M+1=D_j,m for all j∈𝒱  <cit.>.
    Thus this scenario can be represented by an extended graph with one more vertex and an extended 𝐃 with dimension (M+1)×(M+1).




Correspondingly, we define the selection variable 𝐯=[v_1,⋯,v_M]^T∈{0,1}^M, where v_m=1 if the vertex m appears in the path 𝒬 and v_m=0 otherwise.
Furthermore, we define a matrix 𝐖=[W_1,1,⋯,W_1,M;⋯;W_M,1,⋯,W_M,M]∈{0,1}^M× M, with W_y_j,y_j+1=1 for all j=1,⋯,Q-1 and zero otherwise.


With the moving time from the vertex m to the vertex j being D_m,j/a where a is the velocity,
the total moving time along path 𝒬 is

    1/a∑_m=1^M∑_j=1^MW_m,jD_m,j=Tr(𝐃^T𝐖)/a.

Furthermore, since the total motion energy E of the UGV is proportional to the total motion time <cit.>, the motion energy can be expressed in the form of

    E=(α_1/a+α_2)Tr(𝐃^T𝐖),

where α_1 and α_2 are parameters of the model (e.g., for a Pioneer 3DX robot in Fig. 1, α_1=0.29 and α_2=7.4 <cit.>).




 §.§ Backscatter Data Collection Model


Based on the mobility model, the UGV moves along the selected path 𝒬 to collect data from users as shown in Fig. 2.
In particular, from the starting point y_1, the UGV stops for a duration u_y_1 and then it moves along edge (y_1,y_2) to its outward neighbor y_2, and stops for a duration u_y_2. The UGV keeps on moving and stopping along the path until it reaches the destination y_Q.


When the UGV stops at the vertex m (with v_m=1), it will wait for a time duration u_m for data collection.
Out of this u_m, a duration of t_k,m will be assigned to collect data from user k via full-duplex backscatter communication[When user k adapts the variable impedance for modulating the backscattered waveform with information bits, other users keep silent to avoid collision <cit.>.] <cit.>. More specifically, if t_k,m=0, the IoT user k will not be served in duration u_m.
On the other hand, if t_k,m≠0, the RF source at the UGV transmits a symbol x_k,m∈ℂ with 𝔼[|x_k,m|^2]=p_k,m, where p_k,m is the transmit power of the RF source.
Then the received signal-to-noise ratio (SNR) at the UGV tag reader is
η |g_k,m|^2|h_k,m|^2 p_k,m/N_0, where h_k,m∈ℂ is the downlink channel from the UGV to user k, g_k,m∈ℂ is the uplink channel from user k to the UGV, and N_0 is the power of complex Gaussian noise (including the self-interference due to full-duplex backscatter <cit.>).
Furthermore, η is the tag scattering efficiency determined by the load impedance Z_L and the antenna impedance Z_A <cit.>.
For example, in the on-off keying backscatter shown in Fig. 2, the IoT device switches between two load impedances Z_1 and Z_2 with Z_1≠ Z_A and Z_2=Z_A.
This means that the IoT device transmits (Z_L-Z_A)/(Z_L+Z_A)=(Z_1-Z_A)/(Z_1+Z_A) when switching to Z_L=Z_1 and transmits (Z_L-Z_A)/(Z_L+Z_A)=0 when switching to Z_L=Z_2.
Therefore, η=|(Z_1-Z_A)/(Z_1+Z_A)|^2 in the on-off keying.


Based on the backscatter model, the transmission rate during t_k,m is given by

    R_k,m=log_2(1+v_m·βη|g_k,m|^2|h_k,m|^2 p_k,m/N_0),

where β is the performance loss due to imperfect modulation and coding schemes in backscatter communication <cit.>.
For example, in bistatic backscatter communication with frequency shift keying, β=0.5 <cit.>.
On the other hand, in ambient backscatter communication with on-off keying, β is obtained by fitting log_2(1+β x) to 1-ℚ(√(x)) <cit.>, where
ℚ(x)=1/√(2π)∫_x^∞exp(-u^2/2 )du refers to the Q-function.




Remark 1:
The channels {g_k,m,h_k,m} can be pre-determined as follows.
If the environment is static, all the information about object positions, geometry and dielectric properties in the environment is available.
In such a case, ray tracing methods <cit.> could be used to estimate {g_k,m,h_k,m}.
On the other hand, if the channel is varying but with a fixed distribution, we could allow the UGV to collect a small number of measurements at the stopping points before a set of new missions (e.g., five missions) <cit.>.
Then based on the probabilistic framework in <cit.>, the UGV can predict the channels at all the stopping points.



Remark 2:
In practice, the backscatter efficiency would vary with the incident power.
However, if the change of backscatter efficiency is not significant, it is possible to adopt a constant η to facilitate the analysis <cit.>.
On the other hand, if η is not a constant, according to <cit.>, it is possible to compute the range of η.
Then, the worst case approach, which replaces η in (<ref>) with its lower bound η_lb, can be adopted to guarantee the data collection targets for all users.
Notice that this worst case approach can only achieve suboptimal performance, and modeling η as a nonlinear function of the incident power is an important future work.




§ JOINT MOBILITY MANAGEMENT AND POWER ALLOCATION



In wireless data collection systems, the task is to collect certain amount of data from different IoT devices by planning the path (involving variables 𝐯 and 𝐖) and designing the stopping time {t_k,m} and transmit power {p_k,m}.
In particular, the data collection QoS requirement of the k^th IoT device can be described by

    ∑_m=1^M
    t_k,m·log_2(1+v_m·βη|g_k,m|^2|h_k,m|^2 p_k,m/N_0)≥γ_k,

where γ_k>0 (in bit/Hz) is the amount of data to be collected from user k.


Notice that the variables 𝐯 and 𝐖 are dependent since v_m=0 implies W_m,j=W_j,m=0 for any j∈𝒱.
On the other hand, the UGV would visit the vertex with v_m=1, making ∑_j=1^M W_m,j = ∑_j=1^M W_j,m=1.
Combining the above two cases, we have

    ∑_j=1^MW_m,j=v_m, ∑_j=1^MW_j,m=v_m,  ∀ m=1,⋯,M.

Furthermore, since the path must be connected, the following subtour elimination constraints are required to eliminate disjointed sub-tours <cit.>:

    λ_m-λ_j+(∑_l=1^Mv_l-1)W_m,j+(∑_l=1^Mv_l-3)W_j,m
       ≤∑_l=1^Mv_l-2
    +J(2-v_m-v_j)
    ,  ∀ m,j≥ 2, m≠ j,
    
       v_m≤λ_m≤(∑_l=1^Mv_l-1)v_m,  ∀ m≥ 2,

where {λ_m} are slack variables to guarantee a connected path, and ∑_l=1^Mv_l is the number of vertices involved in the path.
The constant J=10^6 is large enough such that the first line of constraint is always satisfied when v_m=0 or v_j=0.
In this way, the vertices not to be visited would not participate in subtour elimination constraints.


Having the data collection and graph mobility constraints satisfied, it is then crucial to reduce the total energy consumption at the UGV.
As the energy consumption includes motion energy
(α_1/a+α_2)Tr(𝐃^T𝐖)
and communication energy ∑_m=1^M∑_k=1^Kt_k,mp_k,m, the joint mobility management and power allocation problem of the data collection system is formulated as[The considered system adopts semi-passive backscatter communication, where the local circuits at IoT devices are powered by their own batteries <cit.>. On the other hand, if the users also request data from the UGV, simultaneous wireless information and power transfer from the UGV to IoT users can be adopted <cit.>.]:


    P1:
       min_𝐯,𝐖,{λ_m}
    {t_k,m,p_k,m}  μ(α_1/a+α_2)Tr(𝐃^T𝐖)
    
                     
    +(2-μ)∑_m=1^M∑_k=1^Kt_k,mp_k,m
    s.t.    ∑_m=1^M
    t_k,m·log_2(1+v_m·βη|g_k,m|^2|h_k,m|^2 p_k,m/N_0)
    
       ≥γ_k,  ∀ k, 
       1/aTr(𝐃^T𝐖)+∑_m=1^M∑_k=1^Kt_k,m≤ T, 
       ∑_j=1^MW_m,j=v_m, ∑_j=1^MW_j,m=v_m,  ∀ m, 
       λ_m-λ_j+(∑_l=1^Mv_l-1)W_m,j+(∑_l=1^Mv_l-3)W_j,m
       ≤∑_l=1^Mv_l-2+J(2-v_m-v_j),
    
       ∀ m,j≥ 2, m≠ j, 
       v_m≤λ_m≤(∑_l=1^Mv_l-1)v_m,  ∀ m≥ 2, 
       W_m,j∈{0,1},  ∀ m,j,  W_m,m=0,  ∀ m, 
       v_1=1, v_m∈{0,1},  ∀ m≥ 2, 
       (1-v_m)· t_k,m=0,  ∀ k,m, 
       t_k,m≥0, p_k,m≥0,  ∀ k,m,


where (<ref>) is for constraining the operation (including moving and data collection) to be completed within T seconds, and (<ref>) is for constraining the stopping time to be zero if the vertex is not visited.
Notice that 0<μ≤ 1 is a weighting factor to control the relative importance between motion energy and communication energy.
Nominally, if we are only interested in minimizing the total energy, we can set μ=1.
On the other hand, if we want to restrict the interference to other co-existing wireless systems, we might set μ<1.


It can be seen from the constraint (<ref>) of P1 that the UGV can choose the stopping vertices, which in turn affect the channel gains to and from the IoT users.
By choosing the stopping vertices with better channel gains to IoT users, the transmit powers {p_k,m} might be reduced.
However, this might also lead to additional motion energy, which in turn costs more energy consumption at the UGV.
Therefore, there exists a trade-off between moving and communication, and solving P1 can concisely balance this energy trade-off.

Unfortunately, problem P1 is nontrivial to solve due to the following reasons.
Firstly, it is NP-hard, since it involves the integer constraints (<ref>)-(<ref>) <cit.>.
Secondly, the data-rate and the energy cost at each vertex are dependent on the transmit power {p_k,m} and transmission time {t_k,m}, which are unknown (see Table I).
This is in contrast to traditional integer programming problems <cit.>, where the reward of visiting each vertex is a constant.






§ OPTIMAL SOLUTION TO P1


Despite the optimization challenges, this section proposes an algorithm that theoretically obtains the optimal solution to P1.
The idea of this algorithm is to eliminate the variables 𝐖,{λ_m},{t_k,m,p_k,m} so as to transform P1 into an equivalent problem only related to 𝐯.
By doing so, we can capitalize on the branch and bound (B&B) method and obtain the optimal solution by pruning out impossible candidates.
In the following, the optimality condition of P1 will be first discussed, which helps in reducing the dimension of P1.




 §.§ Optimality Condition


To address the challenges for solving P1, we first establish the optimality condition of P1.
In particular, by defining

    A_k,m=βη|g_k,m|^2|h_k,m|^2/N_0,

the following proposition (proved in Appendix A) can be established.

The optimal {𝐯^*,t^*_k,m,p_k,m^*} to P1 satisfies:

(i) If m=argmax_l∈𝒱 v_l^*A_k,l, then t^*_k,m≠0; otherwise t^*_k,m=0.

(ii) If t_k,m^*≠ 0, then p^*_k,m≠0.



Proposition 1 indicates that the UGV only needs to allocate time to user k at a single vertex, which is given by m=argmax_l∈𝒱 v^*_lA_k,l.
For other vertices, the allocated time to user k should be zero.
Based on part (i) of Proposition 1, we can set the transmit time

    t_k,m={   s_k
    ,   if m=argmax_l∈𝒱 v_lA_k,l   
       0,   if m≠argmax_l∈𝒱 v_lA_k,l   .,

where s_k>0, without changing the optimal solution to P1.
Correspondingly, the transmit power {p_k,m} can be set to[If m≠argmax_l∈𝒱 v_lA_k,l, we have t_k,m=0. Thus t_k,mp_k,m=0 in the objective of P1 and
t_k,mlog_2(1+v_mA_k,mp_k,m)=0 in (<ref>), meaning that p_k,m would not participate in problem P1.
As a result, we can set p_k,m=0.]

    p_k,m={   q_k
    ,   if m=argmax_l∈𝒱 v_lA_k,l   
       0,   if m≠argmax_l∈𝒱 v_lA_k,l   .,

where q_k is the transmit power corresponding to s_k.
Since s_k>0, by part (ii) of Proposition 1, we also have q_k>0.
Putting (<ref>) and (<ref>) into P1, problem P1 is transformed into


    P2:
    min_𝐯,𝐖,{λ_m}
    {s_k>0,q_k>0}     μ(α_1/a+α_2)Tr(𝐃^T𝐖)
    
       
    +(2-μ)∑_k=1^Ks_kq_k
       s.t.         s_k·log_2[1+(max_l∈𝒱 v_lA_k,l)q_k]
    
       ≥γ_k,  ∀ k, 
       Tr(𝐃^T𝐖)/a+∑_k=1^Ks_k≤ T, 
       (<ref>)-(<ref>).


Notice that  the constraint (<ref>) is dropped since (<ref>) is always satisfied when {t_k,m} takes the form of (<ref>).


The problem P2 is still nontrivial to solve due to the nonlinear coupling between 𝐯 and {𝐖,s_k,q_k} as observed from the constraints (<ref>) and (<ref>).
To resolve such coupling, a straightforward idea is to use alternating minimization for optimizing 𝐯, 𝐖 and {s_k,q_k} iteratively.
However, due to the discrete nature of 𝐯 and 𝐖, such a method could fail to converge.
To this end, this paper proposes to simplify the problem based on elimination of variables.
In particular, we will first derive the optimal solution of 𝐖 and {s_k,q_k} to P2 with fixed 𝐯.
By representing 𝐖 and {s_k,q_k} as functions of 𝐯, problem P2 is simplified to an equivalent problem only involving 𝐯.
Then we will step further to find the optimal solution of vertex selection variable 𝐯.






 §.§ Optimal Solution of 𝐖 and {s_k,q_k} with Fixed 𝐯


When 𝐯=𝐯, where 𝐯 is any feasible solution to P2, the constraint (<ref>) can be dropped since it only involves 𝐯.
On the other hand, the term max_l∈𝒱 v_l A_k,l in constraint (<ref>) becomes a constant, and we denote it as B_k(𝐯):=max_l v_l A_k,l.
Furthermore, it can be seen from the objective function of P2 that q_k is a variable to be minimized and s_klog_2[1+B_k(𝐯)q_k] is a strictly increasing function of q_k.
As a result, the optimal solution of q_k must activate the constraint (<ref>) of P2, which leads to

    q_k=1/B_k(𝐯)(2^γ_k/s_k-1).

Putting (<ref>) into P2, it is proved in Appendix B that the optimal 𝐖^* and {s_k^*} to P2 must activate the constraint (<ref>), i.e., Tr(𝐃^T𝐖^*)/a+∑_k=1^Ks_k^*=T.
Using this result, the quantity Tr(𝐃^T𝐖) in the objective function of P2 can be replaced by a(T-∑_k=1^Ks_k) without changing the problem, and the objective function would be independent of 𝐖.
Then problem P2 is equivalently transformed into the following two-stage optimization problem (detailed procedure given in Appendix B):

    P3:
    min_{s_k>0}    μ(α_1+α_2a)(T-∑_k=1^Ks_k)
    
       
    +(2-μ)∑_k=1^Kγ_k/B_k(𝐯)Θ(s_k/γ_k)
    
      s.t.      ∑_k=1^Ks_k
    =max_𝐖,{λ_m}{T-Tr(𝐃^T𝐖)/a:
    
                   (<ref>)-(<ref>)},

where Θ(x):=x(2^1/x-1).

To solve P3, we first need to compute the right hand side of the constraint, which leads to the following problem:

    max_𝐖,{λ_m} T-Tr(𝐃^T𝐖)/a,     s.t. (<ref>)-(<ref>).

The problem (<ref>) is a travelling salesman problem, which can be optimally solved by the one-tree relaxation algorithm via the software Mosek <cit.>.
In particular, the one-tree relaxation algorithm is an iterative procedure that finds a sequence of one-tree upper bounds to the problem (<ref>) until convergence <cit.>, and the converged solution is guaranteed to be optimal.

Denoting the optimal solution to the problem (<ref>) as {𝐖,λ_m},
the optimal objective value of the travelling salesman problem is given by

    Υ(𝐯):=
    T-Tr(𝐃^T𝐖)/a.

Now, by putting the obtained Υ(𝐯) into P3, the constraint of P3 is written as ∑_k=1^Ks_k
= Υ(𝐯).
Assigning a Lagrange multiplier ρ to this constraint, the Lagrangian of P3 is

    ℒ({s_k},ρ)=   μ(α_1+α_2a)[T-Υ(𝐯)]
    
       
    +(2-μ)∑_k=1^Kγ_k/B_k(𝐯)Θ(s_k/γ_k)
    
       +ρ(∑_k=1^Ks_k-Υ(𝐯)).

According to the first-order KKT condition ∂ℒ/∂ s_k|_s_k=s_k=0, the optimal {s_k} and ρ should together satisfy <cit.>:

    -∇Θ(s_k/γ_k)=
    B_k(𝐯)ρ/2-μ,

where -∇Θ(x) is the gradient of -Θ(x) and is given by

    -∇Θ(x)=1+ln2· 2^1/x/x-2^1/x,

with x>0.
Moreover, as -∇^2Θ(x)=-
ln^22· 2^1/x/x^3<0 and
-∇^3Θ(x)=ln^32· 2^1/x/x^5+3ln^22· 2^1/x/x^4>0, it can be seen that -∇Θ(x) is a strictly decreasing and convex function of x.
As a result, there must exist a strictly decreasing and convex function Λ(x) such that Λ(-∇Θ(x))=x.
That is, the function Λ(x) is the inverse function of -∇Θ(x), and it can be numerically computed and stored as a look-up table, with its shape shown in Fig. 3.
Applying Λ(x) to both sides of (<ref>), we have

    s_k=γ_k·Λ(B_k(𝐯)ρ/2-μ).





Notice that in (<ref>), the only unknown is ρ, which should satisfy the equality constraint of P3:

    ∑_k=1^Kγ_k·Λ(B_k(𝐯)ρ/2-μ)
    =Υ(𝐯).

Since Λ(x) is a decreasing function, bisection search method can be used to find ρ efficiently.
In order to determine the bisection interval, the following proposition (proved in Appendix C) can be established.

The quantity ρ in (<ref>) is bounded as

    -
    (2-μ)∑_k=1^Kγ_k/∑_k=1^Kγ_kB_k(𝐯)∇Θ(Υ(𝐯)/∑_k=1^Kγ_k)
    
                      ≤ρ≤
    -2-μ/min_l B_l(𝐯) ∇Θ(Υ(𝐯)/∑_k=1^Kγ_k).


Once ρ is obtained, we can put ρ into (<ref>) to get {s_k}.
Further putting {s_k} into (<ref>), we have
q_k=(2^γ_k/s_k-1)
/B_k(𝐯).





 §.§ Optimal Solution of 𝐯


With path 𝐖, transmit times {s_k}, and transmit powers {q_k} derived in Section IV-B, the optimal objective value of P3 (equivalently P2 with 𝐯=𝐯) is given by

    Ξ(𝐯)=   μ(α_1+α_2a)[T-Υ(𝐯)]
    
       
    +(2-μ)∑_k=1^Kγ_k/B_k(𝐯)Θ[Λ(B_k(𝐯)
    ρ/2-μ)].

Therefore, problem P2 is re-written as

    P4:min_𝐯     Ξ(𝐯)
    
      s.t.     v_1=1, v_m∈{0,1},  ∀ m=2,⋯,M.

To solve P4, a naive way is to apply exhaustive search for 𝐯.
Unfortunately, since the searching space of {v_m} is very large (i.e., 2^M-1), direct implementation of exhaustive search is impossible.
To address the above issue, a B&B method is presented for systematically pruning out impossible solutions of 𝐯, leading to significant reduction of the computational complexity compared to exhaustive search while guaranteeing the global optimality <cit.>.


In particular, we define the living pool as a set 𝒴 which stores all the solutions that have not been explored, and the incumbent I as the current best objective value that has been obtained.
Initially, I can be set to the objective value of any feasible solution.
On the other hand, since the feasible set of 𝐯 for P4 is Ω={𝐯∈{0,1}^M: v_1=1},
the initial living pool can be set to 𝒴={ℱ_(1,0),ℱ_(1,1)} (notice that 𝒴 is a family of sets over Ω), where
ℱ_(1,0)={𝐯∈Ω: v_2=0} and ℱ_(1,1)={𝐯∈Ω: v_2=1} as shown in Fig. 4a.
It can be seen that Ω=ℱ_(1,0)⋃ℱ_(1,1) and ℱ_(1,0),ℱ_(1,1)⊂Ω.




At the beginning, the B&B method computes a lower bound Ψ(ℱ_(1,0)) such that Ψ(ℱ_(1,0))≤Ξ(𝐯) for any 𝐯∈ℱ_(1,0), where Ξ(𝐯) is defined in (<ref>).
Based on the bounding function value Ψ(ℱ_(1,0)), we consider the following three cases.


  (i) If Ψ(ℱ_(1,0))>I, then the subset ℱ_(1,0) can be discarded, since no feasible solution inside ℱ_(1,0) leads to better objective than the incumbent. Therefore, we update 𝒴←𝒴∖ℱ_(1,0). This case is shown in Fig. 4b.



  (ii) If Ψ(ℱ_(1,0))≤ I and |ℱ_(1,0)|>1, the possibility of a better solution in ℱ_(1,0) cannot be ruled out. As a result, we need to branch on ℱ_(1,0) and generates two subsets
ℱ_(1,0,0)={𝐯∈ℱ_(1,0): v_3=0} and ℱ_(1,0,1)={𝐯∈ℱ_(1,0): v_3=1}.
By treating ℱ_(1,0,0) and ℱ_(1,0,1) as new subset nodes, we update 𝒴←𝒴∖ℱ_(1,0)∪{ℱ_(1,0,0),ℱ_(1,0,1)}.
This case is shown in Fig. 4c.


  (iii) Otherwise, we must have |ℱ_(1,0)|=1. Denoting the unique element in |ℱ_(1,0)| as 𝐯,
we compute Ξ(𝐯).
If Ξ(𝐯)>I, ℱ_(1,0) is discarded.
On the other hand, if Ξ(𝐯)≤ I, the incumbent is updated as I←Ξ(𝐯) and the current best solution is updated as 𝐯^♢←𝐯.
In both cases, since there is no more element in ℱ_(1,0) to be evaluated, we can update 𝒴←𝒴∖ℱ_(1,0).


After ℱ_(1,0) is evaluated, we then evaluate ℱ_(1,1),ℱ_(1,0,0),ℱ_(1,0,1),⋯ by repeating the above procedure until the living pool 𝒴 becomes empty.

For the above B&B method, the key step is to derive the bounding function Ψ.
In particular, to guarantee that the proposed B&B method finds the optimal solution of 𝐯 to P4,
the bounding function needs to satisfy <cit.>:

    Ψ(ℱ_(z_1,⋯,z_N))≤Ξ(𝐯),  ∀ 𝐯∈ℱ_(z_1,⋯,z_N),

where (z_1,⋯,z_N) are the values assigned to (v_1,⋯,v_N) and N is the number of fixed elements in 𝐯.
To this end, consider the following bounding function Ψ:

    Ψ(ℱ_(z_1,⋯,z_N))=   μ(α_1/a+α_2)[T-Φ(𝐳)]
    
       
    +(2-μ)∑_k=1^Kγ_k/B_k([𝐳^T,1_M-N^T]^T)
       ×Θ[Λ(B_k([𝐳^T,1_M-N^T]^T)
    δ/2-μ)],

where 𝐳=[z_1,⋯,z_N]^T∈{0,1}^N (corresponding to the sequence (z_1,⋯,z_N)), and δ is the solution to

    ∑_k=1^K
    γ_k·Λ(B_k([𝐳^T,1_M-N^T]^T)
    δ/2-μ)=Φ(𝐳).

The function

    Φ(𝐳)=   max_𝐖 {T-Tr(𝐃^T𝐖)/a:
    
       ∑_j=1^MW_m,j=∑_j=1^MW_j,m=z_m,
      ∀ m=1,⋯,N, 
       W_m,j∈{0,1},  ∀ m,j,  W_m,m=0,  ∀ m}

represents a bipartite matching problem, which can be numerically computed via the Hungarian algorithm <cit.>.
It is proved in Appendix D that (<ref>) satisfies the property (<ref>).
As a result, by applying the bounding function Ψ in (<ref>), the proposed B&B method is guaranteed to obtain the optimal solution of 𝐯 to problem P4 (equivalently P2).






 §.§ Summary of Algorithm and Complexity Analysis


Since the B&B algorithm finds the optimal solution of 𝐯 to P2, and the optimal solution of 𝐖 and {s_k,q_k} with fixed 𝐯 can be computed according to Section IV-B, the entire algorithm for computing the optimal solution to P2 (equivalently P1) is summarized in Algorithm 1.


In terms of computational complexity, computing Φ(𝐳) in (<ref>) via the Hungarian algorithm requires a complexity of O((2M)^3) <cit.>.
On the other hand, computing δ would require bisection search in solving (<ref>) and the number of iterations is given by log_2(C/ϵ) <cit.>, where C is the length of the initial searching interval given by (<ref>) and ϵ is the target accuracy.
Therefore, computing the bounding function Ψ requires a complexity of O((2M)^3+Klog_2(C/ϵ)).
Finally, computing Ξ(𝐯) for a fixed 𝐯=𝐯 would involve the travelling salesman problem, which requires a complexity of O((M-1)^2·2^M-1) in the worst case <cit.>.

Based on the above analysis, the total complexity of Algorithm 1 is given by

    Comp=   O[X_1((2M)^3
    +Klog_2(C/ϵ))
    
       
    +X_2(M-1)^2·2^M-1],

where X_1 is the number of times for computing Ψ and X_2 is the number of times for computing Ξ(𝐯).
In contrast, the computational complexity for exhaustive search of 𝐯 is 2^M-1·(M-1)^2·2^M-1.
As (2M)^3+Klog_2(C/ϵ) is much smaller than (M-1)^2· 2^M-1, and by simulation X_1+X_2 is significantly smaller than 2^M-1, the proposed Algorithm 1 can significantly reduce the computational complexity compared to exhaustive search.

Finally, for the above algorithm, we need an initial incumbent of I.
Theoretically, I can be set to the function value of any feasible solution to P1.
However, since we are minimizing the objective function of P1, a smaller initial incumbent I could help in reducing the size of the living pool <cit.>, and the next section will derive an efficient initialization method.

Remark 3: If the transmit power is limited due to hardware constraints (e.g., cost of an amplifier), we can add a constraint s_k≤ P_max (with P_max being the upper limit of transmit power) for all k to problem P3.
In such a case, the proposed Algorithm 1 is still applicable if we modify the objective function Ξ(𝐯) in problem P4 into

    Ξ(𝐯)=   min_{s_k} {μ(α_1+α_2a)(T-Υ(𝐯))
    
       
    +(2-μ)∑_k=1^Kγ_k/B_k(𝐯)Θ(s_k/γ_k):
    
       ∑_k=1^Ks_k
    =Υ(𝐯),  0<s_k≤ P_max,  ∀ k
    },

and the bounding function Ψ in (<ref>) into

    Ψ(ℱ_(z_1,⋯,z_N))=min_{s_k} {μ(α_1+α_2a)(T-Φ(𝐳))
    
       
    +(2-μ)∑_k=1^Kγ_k/B_k([𝐳^T,1_M-N^T]^T)Θ(s_k/γ_k):
    
       ∑_k=1^Ks_k
    =Φ(𝐳),  0<s_k≤ P_max,  ∀ k
    }.





§ INITIALIZATION VIA SUCCESSIVE LOCAL SEARCH


In order to obtain a good initial incumbent of I, this section proposes a local optimal solution method based on successive local search <cit.>.
More specifically, we start from a feasible solution of 𝐯 (e.g., 𝐯^[0]=[1,0,⋯,0]^T), and randomly selects a candidate solution 𝐯' from the neighborhood 𝒩(𝐯^[0]). Since a natural neighborhood operator for binary optimization is to flip the value of {v_m}, 𝒩(𝐯^[0]) can be set to

    𝒩(𝐯^[0])={𝐯:||𝐯-𝐯^[0]||_0≤ L, 𝐯∈Ω},

where L≥ 1 is the size of neighborhood.
It can be seen that 𝒩(𝐯^[0]) is a subset of the entire feasible space Ω and containing solutions “close” to 𝐯^[0].


With the neighborhood 𝒩(𝐯^[0]) defined above and the choice of 𝐯 fixed to 𝐯=𝐯', we consider two cases.


  (i) If Ξ(𝐯')≤Ξ(𝐯^[0]), we update 𝐯^[1]←𝐯'. By treating 𝐯^[1] as a new feasible solution, we can construct the next neighborhood 𝒩(𝐯^[1]).


  (ii) If Ξ(𝐯')>Ξ(𝐯^[0]), we find another point within the neighborhood 𝒩(𝐯^[0]) until Ξ(𝐯')≤Ξ(𝐯^[0]).


The above procedure is repeated to generate a sequence of {𝐯^[1],𝐯^[2],⋯} and the converged point is guaranteed to be a local optimal solution to P1 <cit.>.
But since our aim is to obtain a good initial incumbent, it is not necessary to wait until the successive local search converges.
In fact, we can terminate the iterative procedure when the number of iterations is larger than Iter.
Denoting the solution after Iter iterations as 𝐯^♢, we can set the initial incumbent as I=Ξ(𝐯^♢).
The entire procedure to generate an initial incumbent is summarized in Algorithm 2, and the complexity of executing Algorithm 2 is Iter·(M-1)^2·2^M-1.






§ SIMULATION RESULTS AND DISCUSSIONS


This section provides simulation results to evaluate the performance of the UGV backscatter communication system.
It is assumed that the backscattering efficiency is η=0.78 (corresponding to 1.1 dB loss <cit.>), the performance loss due to imperfect modulation is β=0.5 <cit.>, and the weighting factor μ=1.
Within the time budget T=500 s, the data collection targets γ_k∼𝒰(1,2) in the unit of bit/Hz are requested by K=50 IoT users (corresponding to a spectral efficiency of Kγ_k/T=0.1∼0.2 bps/Hz <cit.>), where 𝒰(a,b) is the uniform distribution within the interval [a,b].


Based on the above settings, we simulate the data collection map as a 20 m×20 m=400 m^2 square area, which is a typical size for smart warehouses.
Inside this map, K=50 IoT users[Notice that K is the number of IoT users assigned to the considered UGV. There might exist other IoT users, which can be inactive or assigned to other UGVs.] and M=12 vertices representing stopping points are uniformly scattered.
Among all the vertices, the vertex m=1 is selected as the starting point of the UGV.
With the locations of all the stopping points and the IoT devices, the distances between each IoT device and stopping point can be computed, and the distance-dependent path-loss model ϱ_k,m=ϱ_0·(d_k,m/d_0)^-2.5 is adopted <cit.>, where d_k,m is the distance from user k to the stopping point m, and ϱ_0=10^-3 is the path-loss at d_0=1 m.
Based on the path-loss model, channels g_k,m and h_k,m are generated according to 𝒞𝒩(0,ϱ_k,m).
Each point in the figures is obtained by averaging over 100 simulation runs, with independent channels and realizations of locations of vertices and users in each run.





In order to verify the effectiveness of Algorithm 1 in Section IV, Fig. 5 shows the number of candidate solutions in the living pool 𝒴 versus the number of iterations (represented by Iter in Algorithm 1) when the receiver noise power N_0=-95 dBm  (corresponding to power spectral density -145 dBm/Hz with 100 kHz bandwidth <cit.>).
It can be seen that even with a naive initial incumbent obtained by setting 𝐯=[1,0,⋯,0]^T, the proposed Algorithm 1 still leads to a significantly faster decrease in the number of candidate solutions than the exhaustive search.
Moreover, by using Algorithm 2 to provide an initialization, the number of iterations for Algorithm 1 (together with 20 iterations of Algorithm 2) to reach zero candidate solution can be further decreased.
Notice that in each iteration,  Algorithm 1 requires a smaller computational complexity than that of the exhaustive search according to Section IV-C.
Therefore, the total complexity of Algorithm 1 is significantly reduced compared to the exhaustive search.





Next, we focus on the energy management performance of Algorithm 1.
In particular, the case of K=50 with M=12 is simulated, and the total energy consumption versus the noise power N_0 is shown in Fig. 6a.
For comparison, we also simulate the scheme with no UGV movement (i.e., the optimal solution to P1 with 𝐯=[1,0,⋯,0]^T) and the scheme with full path visiting all vertices (i.e., the optimal solution to P1 with 𝐯=1).
It can be seen that if the noise power is large, by allowing the UGV to visit all the vertices, it is possible to achieve a significantly lower energy consumption compared to the case of no UGV movement.
However, this conclusion does not hold in the small noise power regime, which indicates that moving is not always beneficial.
Fortunately, the proposed Algorithm 1 can automatically determine whether to move and how far to move.
For example, if the noise power is extremely small (e.g., -120 dBm), the UGV could easily collect the data from IoT users at the starting point.
In such a case, the proposed Algorithm 1 would fix the UGV at the starting point.
This can be seen from Fig. 6a at N_0=-120 dBm, in which Algorithm 1 leads to the same performance as the scheme of no UGV movement.
However, if the noise power is increased to a medium value (e.g., -90 dBm), the total energy is reduced by allowing the UGV to move (with the moving path being the red line shown in Fig. 6b).
On the other hand, if the noise power is large (e.g., -60 dBm), the energy for data collection would be high for far-away users.
Therefore, the UGV should spend more motion energy to get closer to IoT users.
This is the black line shown in Fig. 6b.
But no matter which case happens, the proposed algorithm adaptively finds the best trade-off between spending energy on moving versus on communication, and therefore achieves the minimum energy consumption for all the simulated values of N_0 as shown in Fig. 6a.
Notice that the largest transmit power for the proposed Algorithm 1 in Fig. 6 is 30.8 W (occurred at noise power -60 dBm).
Translating this number to the received power density at 1 m gives 2.45 W/m^2, which is within the requirement (<10 W/m^2) set by the IEEE standard C95.1-2005 <cit.>.


The above Fig. 6 has shown that the noise power N_0 can affect the path obtained from Algorithm 1.
In fact, other parameters such as the time budget T and backscatter efficiency η could also impact the path.
To see this, the case of K=50 with M=12 at noise power N_0=-80 dBm is simulated, and the paths for
(T,η)=(500 s,0.78), (T,η)=(20 s,0.78), and (T,η)=(500 s,0.1) are compared in Fig. 7.
It can be seen from Fig. 7 that the path for (T,η)=(20 s,0.78) involves a smaller moving distance than that for (T,η)=(500 s,0.78).
This is because a smaller T would restrict the constraint (<ref>) of P1, which forces the UGV to reduce its motion time and moving distance.
On the other hand, the path for (T,η)=(500 s,0.1) involves a larger moving distance than that for (T,η)=(500 s,0.78), since a smaller η would deteriorate the communication qualities, which forces the UGV to get closer to IoT users.



In order to evaluate the performance of the proposed algorithm under various QoS requirements, the total energy consumption versus the data collection target γ_1=γ_2=⋯=γ_K in bit/Hz at noise power N_0=-90 dBm, is shown in Fig. 8.
It can be seen that under all the simulated values of data collection targets, the proposed Algorithm 1 always achieves the minimum energy consumption.
Moreover, the blue and black curves intersect on the left hand side, while the red and black curves intersect on the right side.
This means that for a very low data collection target, the proposed algorithm results in a non-moving UGV, and for a very high data collection target,
the UGV would visit all the vertices.
Notice that the red curve in Fig. 6 and Fig. 8 is in fact increasing slowly.
The reason behind such a slow change is that the communication energy is negligible compared to the energy required for mobility if the UGV visits all vertices.







To further assess the performance of the proposed Algorithm 1 when μ varies, the case of K=50 with M=12 at N_0=-90 dBm is simulated and the result is shown in Fig. 9.
It can be seen that if μ=0, the proposed Algorithm 1 has the same performance as that of the full path.
This is because μ=0 would make the motion energy disappear in the objective function of P1, and the best strategy is to allow the UGV to visit all vertices.
However, even with a slight increase in μ, the proposed Algorithm 1 would outperform other benchmark schemes.
Finally, it can be seen that if μ increases, the motion energy decreases but the communication energy increases.
Therefore, μ can be used to adjust the relative amount of communication energy versus motion energy.
Notice that the communication energy of UGV in Fig. 9 is in the same order of magnitude as the motion energy of UGV when μ=1.
This is different from UAV communications, where the propulsion energy of UAV is much larger than the communication energy to keep the UAV aloft <cit.>.





Finally, we analyze the impact of the number of vertices M on the energy consumption.
In particular, the case of K=50 with M∈{3,6,9,12} at noise power N_0=-90 dBm is simulated, and the result is shown in Fig. 10.
It can be seen from Fig. 10 that the scheme with no UGV movement is independent of M, since the UGV is fixed at the starting point and would not visit any other vertex.
On the other hand, the performance of the full path visiting all vertices becomes worse when M increases.
This is because the UGV needs to visit more vertices and therefore consumes more motion energy.
In addition, the total energy consumption of the proposed Algorithm 1 decreases when M increases, since a larger M would give the algorithm more freedom to optimize the trajectory.
Lastly, for a fixed value of M, the performance gap between the best case and the worst case could be large as shown in the error bars of Fig. 10.
This is due to the different spatial distributions of users in different map realizations.
If the users are distributed in clusters (e.g., users' locations follow Gaussian mixture distribution), the proposed algorithm could achieve the best performance.
On the other hand, if the users are widely spread out in the map, it would be difficult to collect data from these users, resulting in the worst case performance of the algorithm.
It is worth noting that the locations of vertices are in general independent of the users' locations.
This is because the vertices should be placed where the UGV is able to approach and stop.
However, if the locations of vertices are allowed to be chosen freely, a promising heuristic for setting the locations of vertices is to cluster the K users into M groups and place the vertices at the cluster centers <cit.>.





§ CONCLUSIONS


This paper studied a UGV-based backscatter data collection system, with an integrated graph mobility model and backscatter communication model.
The joint mobility management and power allocation problem was formulated with the aim of energy minimization subject to communication QoS constraints and mobility graph structure constraints.
An algorithm that achieves the optimal solution was derived, and it automatically balances the trade-off between spending energy on moving and on communication.
Simulation results showed that the proposed algorithm could significantly save energy consumption compared to the scheme with no UGV movement and the scheme with a fixed moving path.





§ PROOF OF PROPOSITION 1


This proposition contains two parts, and we will first prove part (ii) by contradiction.



 §.§ Proof of Part (ii)


To prove part (ii) by contradiction, consider an optimal solution {t_k,m^*,p_k,m^*,𝐯^*,𝐖^*,λ_m^*} to P1 with a particular
(i,j) such that t_i,j^*=t≠ 0.
Assume that the corresponding p_i,j^*=0.

Since {t_k,m^*,p_k,m^*,𝐯^*,𝐖^*,λ_m^*} is optimal to P1, it must satisfy (<ref>) of P1, i.e., there must exist some n≠ j such that t_i,n^*=t'≠0 and p_i,n^*=p'≠0.
Now, consider a related problem of P1 by fixing all the variables to their optimal values except for (t_i,j,t_i,n,p_i,j,p_i,n):

    min_t_i,j,t_i,n,p_i,j,p_i,n≥ 0     t_i,jp_i,j+t_i,np_i,n
          s.t.             C_1+t_i,jlog_2(1+v_j^*A_i,jp_i,j)
    
       
    +t_i,nlog_2(1+v_n^*A_i,np_i,n)≥γ_i, 
       C_2+t_i,j+t_i,n≤ T,

where

    C_1=   ∑_(k,m)∉{(i,j),(i,n)}t_k,m^*log_2(1+v_m^*A_k,mp_k,m^*),
    
    
    C_2=   ∑_(k,m)∉{(i,j),(i,n)}t_k,m^*+Tr(𝐃^T𝐖^*)/a.

As {t_k,m^*,p_k,m^*,𝐯^*,𝐖^*} is optimal to P1, it can be seen that

    (t_i,j,t_i,n,p_i,j,p_i,n)=(t,t',0,p') is optimal to (<ref>).

Therefore, (t,t',0,p') should satisfy the constraints of (<ref>), which leads to

    C_1+t'log_2(1+v_n^*A_i,np')≥γ_i,  
    C_2+t+t'≤ T.

Furthermore, by Jensen's inequality, we have

    (t+t')log_2(1+v_n^*A_i,n·t·0+t'p'/t+t')
    
       
    >
    tlog_2(1+v_n^*A_i,n· 0)+
    t'log_2(1+v_n^*A_i,np'),

where the strict inequality is due to t,t',p' ≠ 0.
Adding C_1 to both sides of (<ref>), and combining the first inequality of (<ref>), we have C_1+(t+t')log_2(1+v_n^*A_i,n·t'/t+t'p')>γ_i.
Comparing this result and the second inequality of (<ref>) to the constraints of (<ref>), it can be seen that (t_i,j,t_i,n,p_i,j,p_i,n)=(0,t+t',0,t'/t+t'p'-Δ p) is feasible for (<ref>) under sufficiently small Δ p>0.
Putting (0,t+t',0,t'/t+t'p'-Δ p) and (t,t',0,p') from (<ref>) into the objective function of (<ref>), we obtain t'p'-(t+t')Δ p and t'p', respectively.
Further due to t'p'-(t+t')Δ p<t'p', it is clear that (t,t',0,p') cannot be optimal to (<ref>).
This contradicts to (<ref>).
Therefore, p_i,j^*≠ 0.




 §.§ Proof of the First Part of (i)


The first part of (i) can be proved by following a similar procedure to that of part (ii).
In particular, consider an optimal solution {t_k,m^*,p_k,m^*,𝐯^*,𝐖^*,λ_m^*} to P1.
Assume that there exists some user i such that t_i,j^*=0 at vertex j=argmax_l∈𝒱 v_l^*A_i,l.

Since {t_k,m^*,p_k,m^*,𝐯^*,𝐖^*,λ_k^*} is optimal to P1, it must satisfy (<ref>) of P1, i.e., there must exist some n≠ j such that t_i,n^*=t'≠0 and p_i,n^*=p'≠0.
Furthermore, under {t_k,m^*,p_k,m^*,𝐯^*,𝐖^*,λ_m^*} with t_i,j^*=0, it can be shown that (t_i,j,t_i,n,p_i,j,p_i,n)=(0,t',p_i,j^*,p') is optimal to (<ref>).
As a result, (0,t',p_i,j^*,p') must satisfy the constraints of (<ref>), i.e.,

    C_1+t'log_2(1+v_n^*A_i,np')≥γ_i,  
    C_2+t'≤ T,

which can be re-written as

    C_1+t'log_2(1+v_j^*A_i,j·v_n^*A_i,n/v_j^*A_i,jp')≥γ_i,  C_2+t'≤ T.

Comparing (<ref>) with the constraints of (<ref>), it can be seen that (t_i,j,t_i,n,p_i,j,p_i,n)=(t',0,v_n^*A_i,n/v_j^*A_i,jp',0) is also feasible for (<ref>).
Putting (0,t',p_i,j^*,p') and (t',0,v_n^*A_i,n/v_j^*A_i,jp',0) into the objective function of (<ref>), we obtain p't' and
v_n^*A_i,n/v_j^*A_i,jp't', respectively.
Further due to (0,t',p_i,j^*,p') being optimal to (<ref>), we have p't'≤v_n^*A_i,n/v_j^*A_i,jp't', which leads to
v_j^*A_i,j≤ v_n^*A_i,n.
As A_i,j and A_i,n are not equal almost surely, and v_j^*,v_n^* ∈{0,1}, the equality sign in v_j^*A_i,j≤ v_n^*A_i,n cannot hold.
This gives us v_j^*A_i,j<v_n^*A_i,n, but it contradicts to j=argmax_l∈𝒱 v_l^*A_i,l.
Therefore, t_i,j^*≠ 0 at vertex j=argmax_l∈𝒱 v_l^*A_i,l.




 §.§ Proof of the Second Part of (i)


We will prove the second part of (i) by contradiction.
In particular, assume that there exists some user i such that t^*_i,n≠ 0 at vertex n≠argmax_l∈𝒱 v_l^*A_i,l.
On the other hand, based on the first part of (i) of Proposition 1, we also have t^*_i,j≠ 0 at vertex j=argmax_l∈𝒱 v_l^*A_i,l.
Correspondingly, by part (ii) of Proposition 1, we have p^*_i,n,p^*_i,j≠ 0.

Now, the partial Lagrangian of P1 with respect to {t_k,m,p_k,m} under fixed 𝐯=𝐯^* is

    ℒ({t_k,m,p_k,m},{ζ_k,φ,χ_k,m,θ_k,m,ξ_k,m})
    
    
    =   (2-μ)∑_m=1^M∑_k=1^Kt_k,mp_k,m
       
    +∑_k=1^K
    ζ_k[γ_k-1/T∑_m=1^M
    t_k,mlog_2(1+v_m^*A_k,mp_k,m)]
    
       
    +φ(∑_m=1^M∑_k=1^Kt_k,m+
    1/aTr(𝐃^T𝐖)-T)
    
       
    +∑_m=1^M∑_k=1^Kχ_k,m[(1-v_m^*)t_k,m]
    
       -∑_m=1^M∑_k=1^Kθ_k,mt_k,m-∑_m=1^M∑_k=1^Kξ_k,mp_k,m
    ,

where {ζ_k,φ,χ_k,m,θ_k,m,ξ_k,m} are Lagrange multipliers.
Since P1 is convex in {t_k,m} with {p_k,m} fixed (vice versa), according to the Karush-Kuhn-Tucker condition <cit.>, the optimal {t^*_k,m,p^*_k,m} and {ζ_k^*,φ^*,χ_k,m^*,θ_k,m^*,ξ_k,m^*}
must satisfy


    θ^*_k,mt^*_k,m=0, ξ^*_k,mp^*_k,m=0,  ∀ k,m, 
       
    (2-μ)t_k,m^*-ζ^*_kt^*_k,m/Tln2·v_m^*A_k,m/1+v_m^*A_k,mp_k,m^*-ξ^*_k,m
       
    =0,  ∀ k,m, 
       ζ^*_klog_2(1+v_m^*A_k,mp^*_k,m)-(2-μ)Tp_k,m^*
    
       
    =T[φ^*-θ^*_k,m+χ_k,m^*(1-v_m^*)],  ∀ k,m.


Putting t^*_i,j,t^*_i,n,p_i,j^*,p_i,n^*≠ 0 into (<ref>), we have

    θ^*_i,j=θ^*_i,n=ξ^*_i,j=ξ^*_i,n=0.

Further putting ξ^*_i,n=0 from (<ref>) into (<ref>), the following equation is obtained

    p^*_i,n=ζ_i^*/ln2· (2-μ)T-1/v_n^*A_i,n.

Substituting (<ref>), θ^*_i,n=0 from (<ref>), and v_n^*=1 (due to (<ref>) and t^*_i,n≠ 0) into (<ref>), equation (<ref>) is reformulated as
F(v_n^*A_i,n)=Tφ^*, where

    F(x)=
    ζ^*_ilog_2[ζ_i^*x/ln2· (2-μ)T]
    +(2-μ)T/x
    -ζ_i^*/ln2

with x≠ 0.
Notice that due to v_n^*A_i,n≠ 0 (as v_n^*=1), F(v_n^*A_i,n) is well-defined.
Similarly, by using ξ^*_i,j=θ^*_i,j=0 from (<ref>), we obtain F(v_j^*A_i,j)=Tφ^*. Therefore,

    F(v_n^*A_i,n)=F(v_j^*A_i,j).

Now, the derivative of F(x) can be computed to be

    ∇_xF=   1/x[ζ^*_i/ln2-(2-μ)T/x]
    
    
    =   1/x[(2-μ)Tp^*_i,n+(2-μ)T(1/v_n^*A_i,n-1/x)],

where the second equality is obtained from (<ref>).
Since p^*_i,n>0, it is clear that ∇_xF(x)>0 for any x∈[v^*_nA_i,n,v^*_jA_i,j] and F(x) is a strictly increasing function of x over this interval.
Combining the result from (<ref>), we have v^*_nA_i,n=v_j^*A_i,j=max_l∈𝒱 v^*_lA_i,l.
This contradicts to n≠argmax_l∈𝒱 v_l^*A_i,l.
Therefore, t^*_i,n=0 at vertex n≠argmax_l∈𝒱 v_l^*A_i,l.



§ TRANSFORMATION FROM P2 TO P3 WITH FIXED 𝐯=𝐯


With q_k in (<ref>), the first constraint of P2 is always satisfied.
Therefore, we can re-write P2 as

    min_{s_k>0,𝐖,λ_m}     μ(α_1/a+α_2)Tr(𝐃^T𝐖)
    
       
    +(2-μ)∑_k=1^Ks_k(2^γ_k/s_k-1)/B_k(𝐯)
        s.t.           ∑_k=1^Ks_k≤ T-
    1/aTr(𝐃^T𝐖),
    
       (<ref>)-(<ref>).

Since s_k is not involved in (<ref>)-(<ref>), the partial Lagrangian of the problem (<ref>) with respect to {s_k} is given by

    ℒ({s_k},φ,{θ_k})=
    (2-μ)∑_k=1^Ks_k(2^γ_k/s_k-1)/B_k(𝐯)
       +φ[∑_k=1^Ks_k-T+1/aTr(𝐃^T𝐖)]
    -∑_k=1^Kθ_ks_k,

where φ≥ 0 and θ_k≥ 0 are Lagrange multipliers.
Since problem (<ref>) is convex in {s_k}, we must have ∂ℒ/∂ s_k=0 according to the KKT condition,
and the optimal {s_k} and (φ,θ_k) should together satisfy

    (2-μ)/B_k(𝐯)·(1+ln2· 2^γ_k/s_k·γ_k/s_k-2^γ_k/s_k)
    =φ-θ_k.

Since the derivative of the left hand side of (<ref>) with respect to s_k is
-(2-μ)/B_k(𝐯)·ln^22· 2^γ_k/s_k·γ_k^2/(s_k)^3<0,
the left hand side function in (<ref>) is a strictly decreasing function of s_k.
Furthermore, due to

    lim_s_k→+∞(2-μ)/B_k(𝐯)·(1+ln2· 2^γ_k/s_k·γ_k/s_k-2^γ_k/s_k)
    =0,

it is clear that φ-θ_k>0.
As θ_ks_k=0 from the complementary slackness condition, and due to s_k>0,
we must have θ_k=0, and therefore φ>0.

Finally, from φ>0 and the complementary slackness condition
φ[1/aTr(𝐃^T𝐖)+∑_k=1^Ks_k-T]=0,
the equation Tr(𝐃^T𝐖)=a(T-∑_k=1^Ks_k) holds.
Substituting this result into P2, P2 becomes

    min_{s_k>0,𝐖,λ_m}    μ(α_1+α_2a)(T-∑_k=1^Ks_k)
    
       
    +(2-μ)∑_k=1^Ks_k/B_k(𝐯)(2^γ_k/s_k-1),
    
    s.t.         ∑_k=1^Ks_k=T-
    1/aTr(𝐃^T𝐖),
    
       (<ref>)-(<ref>).

Now the objective function of (<ref>) is independent of 𝐖 and is a decreasing function of s_k (since its derivative can be shown to be negative).
Therefore, the minimum value of the objective function in (<ref>) is obtained when s_k is maximized.
Since ∑_k=1^Ks_k=T-Tr(𝐃^T𝐖)/a, maximizing T-Tr(𝐃^T𝐖)/a helps in enlarging the maximum values of {s_k}. Therefore, (<ref>) is equivalently transformed into P3.





§ PROOF OF PROPOSITION 2


We first prove the left part of equation (<ref>).
In particular, since Λ(x) is a convex function (see Fig. 3), we must have

    Λ(1/∑_k=1^Kγ_k∑_k=1^Kγ_k·B_k(𝐯)ρ/2-μ)
    
       ≤1/∑_k=1^Kγ_k∑_k=1^Kγ_k·Λ(B_k(𝐯)ρ/2-μ)
    =Υ(𝐯)/∑_k=1^Kγ_k,

where “≤” is due to Jensen's inequlity and “=” is due to (<ref>).
By further applying function -∇Θ to both sides of (<ref>), and since -∇Θ is a strictly decreasing function,
 (<ref>) becomes

    1/∑_k=1^Kγ_k∑_k=1^Kγ_kB_k(𝐯)·ρ/2-μ≥-∇Θ(Υ(𝐯)/∑_k=1^Kγ_k),

which immediately leads to the left part of (<ref>).

Next, we prove the right part of equation (<ref>).
More specifically, since B_k(𝐯)≥min_l B_l(𝐯) and ρ≥0, we have

    B_k(𝐯)ρ/2-μ≥min_l B_l(𝐯)ρ/2-μ.

Applying Λ(x) to both sides of (<ref>), and since Λ(x) is a strictly decreasing function,
 (<ref>) becomes

    Λ(B_k(𝐯)ρ/2-μ)
    ≤Λ(min_l B_l(𝐯)ρ/2-μ).

Based on the above result and equation (<ref>), it is clear that

    ∑_k=1^Kγ_k·Λ(min_l B_l(𝐯)ρ/2-μ)
    ≥Υ(𝐯),

which is equivalent to

    min_l B_l(𝐯)ρ/2-μ≤
    -∇Θ(Υ(𝐯)/∑_k=1^Kγ_k).

This leads to the right part of (<ref>), and the proof is completed.



§ PROOF OF LOWER BOUND PROPERTY OF Ψ IN (<REF>)


To begin with, it is noticed that Ξ(𝐯) in (<ref>) is the optimal value of problem P3.
Therefore, to prove Ψ(ℱ_(z_1,⋯,z_N))≤Ξ(𝐯) for any 𝐯∈ℱ_(z_1,⋯,z_N), we consider two relaxations for P3 with 𝐯∈ℱ_(z_1,⋯,z_N).

(i) Based on the definition B_k(𝐯):=max_l v_l A_k,l, it can be seen that
B_k(𝐯_1)≥ B_k(𝐯_2) if 𝐯_1≽𝐯_2, where ≽ means “Pareto dominance”.
Furthermore, since [𝐳^T,1_M-N^T]^T≽𝐯 for any 𝐯∈ℱ_(z_1,⋯,z_N), we have
B_k([𝐳^T,1_M-N^T]^T)≥ B_k(𝐯) for any 𝐯∈ℱ_(z_1,⋯,z_N).
Using this result, the objective function of P3 can be lower bounded by μ(α_1+α_2a)(T-∑_k=1^Ks_k)
+(2-μ)∑_k=1^Kγ_k/B_k([𝐳^T,1_M-N^T]^T)Θ(s_k/γ_k).

(ii) Notice that

    max_𝐖,{λ_m}{T-1/aTr(𝐃^T𝐖):(<ref>)-(<ref>),𝐯∈ℱ_(z_1,⋯,z_N)}
       ≤max_𝐖{T-1/aTr(𝐃^T𝐖):(<ref>),(<ref>),𝐯∈ℱ_(z_1,⋯,z_N)}
       
    =Φ(𝐳),

where the inequality is obtained by dropping constraints (<ref>) and (<ref>).
Using (<ref>), the constraint of P3 can be relaxed into
∑_k=1^Ks_k≤Φ(𝐳).

Based on the above two relaxations, the problem P3 with 𝐯∈ℱ_(z_1,⋯,z_N) is relaxed into

    min_{s_k>0}     μ(α_1+α_2a)(T-∑_k=1^Ks_k)
    
       
    +(2-μ)∑_k=1^Kγ_k/B_k([𝐳^T,1_M-N^T]^T)Θ(s_k/γ_k)
    
      s.t.       ∑_k=1^Ks_k≤Φ(𝐳).

Using the result from (<ref>)-(<ref>), the optimal {s_k} to the problem (<ref>) is given by

    s_k^*=γ_k·Λ(B_k([𝐳^T,1_M-N^T]^T)δ/2-μ),

with δ obtained from (<ref>).
Putting (<ref>) into the objective function of (<ref>), we immediately obtain Ψ in (<ref>), which is obviously a lower bound to the objective function Ξ(𝐯) of P3 for any 𝐯∈ℱ_(z_1,⋯,z_N).


51

iot1 J. Chen, K. Hu, Q. Wang, Y. Sun, Z. Shi, and S. He, “Narrowband Internet of Things: Implementations and applications,” IEEE IoT Journal, vol. 4, no. 6, pp. 2309-2314, Dec. 2017.

wpt M. Xia and S. Aïssa, “On the efficiency of far-field wireless power transfer,” IEEE Trans. Signal Process., vol. 63, no. 11, pp. 2835-2847, Jun. 2015.

back1 V. Liu, A. Parks, V. Talla, S. Gollakota, D. Wetherall, and J. R. Smith, “Ambient backscatter: Wireless communication out of thin air,” in Proc. ACM SIGCOMM, 2013, pp. 39-50.

back2 J. Kimionis, A. Bletsas, and J. N. Sahalos, “Increased range bistatic scatter radio,” IEEE Trans. Commun., vol. 62, no. 3, pp. 1091-1104, Mar. 2014.

back3 C. Boyer and S. Roy, “Backscatter communication and RFID: Coding, energy, and MIMO analysis,” IEEE Trans. Commun., vol. 62, no. 3, pp. 770-785, Mar. 2014.

back4 A. Alma'aitah, H. S. Hassanein, and M. Ibnkahla, “Tag modulation silencing: Design and application in RFID anti-collision protocols,” IEEE Trans. Commun., vol. 62, no. 11, pp. 4068-4079, Nov. 2014

back5 B. Clerckx, Z. B. Zawawi, and K. Huang, “Wirelessly powered backscatter communications: Waveform design and SNR-energy tradeoff,” IEEE Commun. Lett., vol. 21, no. 10, pp. 2234-2237, Oct. 2017.

back6 B. Lyu, H. Guo, Z. Yang, and G. Gui, “Throughput maximization for hybrid backscatter assisted cognitive wireless powered radio networks,” IEEE IoT Journal, vol. 5 , no. 3, pp. 2015-2024, Jun. 2018.

back7 G. Wang, F. Gao, R. Fan, and C. Tellambura, “Ambient backscatter communication systems: Detection and performance analysis,” IEEE Trans. Commun, vol. 64, no. 11, pp. 4836-4846, Nov. 2016.

back8 J. Qian, F. Gao, G. Wang, S. Jin, and H. Zhu, “Noncoherent detections for ambient backscatter systems,” IEEE Trans. Wireless Commun., vol. 16, no. 3, pp. 1412-1422, Mar. 2017.

ugv2 Y. Mei, Y. H. Lu, Y. Hu, and C. Lee, “Deployment of mobile robots with energy and timing constraints,” IEEE Trans. Robotics, vol. 22, no. 3, pp. 507-522, Jun. 2006.

ugv3 G. Wang, M. J. Irwin, P. Berman, H. Fu, and T. F. L. Porta, “Optimizing sensor movement planning for energy efficiency,” in Proc. ISLPED, pp. 215-220, 2005.

ugv4 Y. Yan and Y. Mostofi, “Co-optimization of communication and motion planning of a robotic operation under resource constraints and in fading environments,” IEEE Trans. Wireless Commun., vol. 12, no. 4, pp. 1562-1572, Apr. 2013.

ugv5 Y. Shu, H. Yousefi, P. Cheng, J. Chen, Y. Gu, T. He, and K. G. Shin, “Near-optimal velocity control for mobile charging in wireless
rechargeable sensor networks,” IEEE Trans. Mobile Comput., vol. 15, no. 7, pp. 1699-1713, Jul. 2016.

ugv6 S. Wang, M. Xia, K. Huang, and Y.-C. Wu, “Wirelessly powered two-way communication with nonlinear energy harvesting model: Rate regions under fixed and mobile relay,” IEEE Trans. Wireless Commun., vol. 16, no. 12, pp. 8190 -8204, Dec. 2017.

laporte1 G. Laporte, “The vehicle routing problem: An overview of exact and approximate algorithms,” European Journal of Operational Research, vol. 59, no. 3, pp. 345-358, 1992.

ugv7 M. Ma, Y. Yang, and M. Zhao, “Tour planning for mobile data gathering mechanisms in wireless sensor networks,” IEEE Trans. Veh. Technol., vol. 62, no. 4, pp. 1472-1483, May 2013.

ugv8 M. Zhao, J. Li, and Y. Yang, “A framework of joint mobile energy replenishment and data gathering in wireless rechargeable sensor networks,” IEEE Trans. Mobile Comput., vol. 13, no. 12, pp. 2689-2705, Dec. 2014.

uav1 Y. Zeng and R. Zhang, “Energy-efficient UAV communication with trajectory optimization,” IEEE Trans. Wireless Commun., vol. 16, no. 6, pp. 3747-3760, Jun. 2017.

uav2 Q. Wu, Y. Zeng, and R. Zhang, “Joint trajectory and communication design for multi-UAV enabled wireless networks,” IEEE Trans Wireless Commun., vol. 17, no. 3, pp. 2109-2121, Mar. 2018.

uav3 Y Sun, D. W. K. Ng, D. Xu, L. Dai, and R. Schober, “Resource allocation for solar powered UAV communication systems,” in Proc. IEEE SPAWC'18, Kalamata, Greece, Jun. 2018.

uav4H. Sallouha, M. M. Azari, and S. Pollin, “Energy-constrained UAV trajectory design for ground node localization,” in Proc. IEEE GLOBECOM, Abu Dhabi, UAE, Dec. 2018.

uav5 H. Sallouha, M. M. Azari, A. Chiumento, and S. Pollin, “Aerial anchors positioning for reliable rss-based outdoor localization in urban environments,” IEEE Wireless Commun. Lett., vol. 7, no. 3, pp.  376-379, Jun. 2018.


graph J. A. Bondy and U. Murthy, Graph Theory with Applications. New York: Elsevier, 1976.

laporte2 G. Laporte, “The traveling salesman problem: An overview of exact and approximate algorithms,” European Journal of Operational Research, vol. 59, no. 2, pp. 231-247, 1992.

FD1 V. Liu, V. Talla, and S. Gollakota, “Enabling instantaneous feedback with full-duplex backscatter,” in Proc. ACM MobiCom, 2014, pp. 67-78.

FD2 W. Liu, K. Huang, X. Zhou, and S. Durrani, “Full-duplex backscatter interference networks based on time-hopping spread spectrum,” IEEE Trans. Wireless Commun., vol. 16, no. 7, pp. 4361-4377, Jul. 2017.

back10 G. Zhu, S. W. Ko, and K. Huang, “Inference from randomized transmissions by many backscatter sensors,” IEEE Trans. Wireless Commun., vol. 17, no. 5, pp. 3111-3127, May 2018.

back11 S. H. Kim and D. I. Kim, “Hybrid backscatter communication for wireless-powered heterogeneous networks,” IEEE Trans. Wireless Commun., vol. 16, no. 10, pp. 6557-6570, Oct. 2017.

ook J. G. Proakis, Digital Communications (4th edition). New York, NY, USA: McGraw-Hill, 2001.

E1  K. A. Remley, H. R. Anderson, and A. Weisshar, “Improving the accuracy of ray-tracing techniques for indoor propagation modeling,” IEEE Trans. Veh. Technol., vol. 49, no. 6, pp. 2350-2358, Nov. 2000.

E3 M. Malmirchegini and Y. Mostofi, “On the spatial predictability of communication channels,” IEEE Trans. Wireless Commun., vol. 11, no. 3, pp. 964-978, Mar. 2012.

back13 P. N. Alevizos, K. Tountas, and A. Bletsas, “Multistatic scatter radio sensor networks for extended coverage,” IEEE Trans. Wireless Commun., vol. 17, no. 7, pp. 4522-4535, Jul. 2018.

back12 A. Bletsas, A. G. Dimitriou, and J. N. Sahalos, “Improving Backscatter Radio Tag Efficiency,” IEEE Trans. Microw. Theory Techn., vol. 58, no. 6, pp. 1502-1509, Jun. 2010.

Claessens1 S. Claessens, D. Schreurs, and S. Pollin, “SWIPT with biased ASK modulation and dual-purpose hardware," in Proc. WPTC, Taipei, May 2017.

Claessens2 S. Claessens, N. Pan, M. Rajabi, D. Schreurs, and S. Pollin, “Enhanced biased ASK modulation performance for SWIPT with AWGN channel and dual-purpose hardware,” IEEE Trans. Microw. Theory Techn., vol. 66, no. 7, pp. 3478-3486, Jul. 2018.

Claessens3 S. Claessens, M. Rajabi, N. Pan, D. Schreurs, and S. Pollin, “Two-tone FSK modulation for SWIPT”, in Proc. WPTC, Montreal, 2018.

Claessens4 M. Rajabi, N. Pan, S. Claessens, S. Pollin, and D. Schreurs, “Modulation techniques for simultaneous wireless information and power transfer With an integrated rectifier receiver,” IEEE Trans. Microw. Theory Techn., vol. 66, no. 5, pp. 2373-2385, May 2018.

Claessens5 D. I. Kim, J. H. Moon, and J. J. Park, “New SWIPT Using PAPR: How it works,” IEEE Wireless Commun. Lett., vol. 5, no. 6, pp. 672-675, Dec. 2016.


dimitri D. P. Bertsekas, Network Optimization: Continuous and Discrete Models. Athena Scientific, 1998.

opt1 S. Boyd and L. Vandenberghe, Convex Optimization. Cambridge, U.K.: Cambridge Univ. Press, 2004.

bb J. Clausen, Branch and Bound Algorithms: Principles and Examples. Copenhagen, Denmark: Univ. Copenhagen, 1999.

MINLP P. Belotti, C. Kirches, S. Leyffer, J. Linderoth, J. Luedtke, and A. Mahajan, “Mixed-integer nonlinear optimization,” Acta Numerica, vol. 22, pp. 1-131, 2013.

hungarian H. W. Kuhn, “The Hungarian method for the assignment problem,” Naval Research Logistics, vol. 2, no. 1-2, pp. 83-97, 1955.

bisection K. E. Atkinson, An Introduction to Numerical Analysis (2nd edition). New York: John Wiley and Sons, 1989.

tsp M. Held and R. M. Karp, “A dynamic programming approach to sequencing problems,” J. of Soc. for Indust. and Appl. Math., vol. 10, no. 1, pp. 196-210, Mar. 1962.

hybrid E. G. Talbi, “Combining metaheuristics with mathematical programming, constraint programming and machine learning,” Annals of Operations Research, vol. 240, no. 1, pp. 171-215, May 2016.

meta M. Gendreau and JY Potvin, Handbook of Metaheuristics (2nd edition). New York: Springer; 2010.

rls F. Neumann and I. Wegener, “Randomized local search, evolutionary algorithms, and the minimum spanning tree problem,” Theoretical Computer Science, vol. 378, no. 1, pp. 32-40, 2007.

channel S. Wang, M. Xia, and Y.-C. Wu, “Multicast wirelessly powered network with large number of antennas via first-order method,” IEEE Trans. Wireless Commun., vol. 17, no. 6, pp. 3781-3793, Jun. 2018.

cluster Y. Yan and Y. Mostofi, “Efficient clustering and path planning strategies for robotic data collection using space-filling curves,” IEEE Transactions on Control of Network Systems, vol. 4, no. 4, pp. 838-849, Dec. 2017.






