Backscatter Data Collection with Unmanned Ground Vehicle: Mobility Management and Power Allocation

TL;DR

This paper proposes an energy-efficient data collection scheme using unmanned ground vehicles with backscatter communication, balancing movement and communication energy to optimize data collection from IoT devices.

Contribution

It introduces a joint mobility and power allocation model for UGV-based backscatter data collection, along with an algorithm achieving minimum energy consumption.

Findings

Lower noise power favors minimal UGV movement.

Higher noise power encourages UGV to move closer to IoT devices.

The proposed algorithm balances movement and communication energy effectively.

Abstract

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.

Tables

Table 1. TABLE I: Summary of Symbol Notations

Variable	Description
$v_m\in \{0,1\}$	$v_m=1$ represents the vertex $m$ being involved in the path; $v_m=0$ otherwise.
$W_{m,j}\in \{0,1\}$	$W_{m,j}=1$ represents the edge $(m,j)$ being involved in the path; $W_{m,j}=0$ otherwise.
$t_{k,m}\in {ℝ}_{+}$	Time (in $s$) allocated to user $k$ when UGV is at the $m^{\mathrm{th}}$ stopping point.
$p_{k,m}\in {ℝ}_{+}$	Transmit power (in $\mathrm{Watt}$) to user $k$ when UGV is at the $m^{\mathrm{th}}$ stopping point.
Parameter	Description
$M$	Number of stopping points.
$K$	Number of users.
$𝒱,ℰ$	$𝒱$ ($ℰ$) is the set of all vertices (edges).
$D_{m,j}$	Distance (in $m$) from the $m^{\mathrm{th}}$ vertex to the $j^{\mathrm{th}}$ vertex.
$\mu$	Weighting factor of motion energy.
${\alpha }_1,{\alpha }_2$	Parameters of the UGV motion energy model.
$a$	Constant velocity (in $m/s$) of the UGV.
$g_{k,m},h_{k,m}$	$g_{k,m}$ ($h_{k,m}$) is the downlink (uplink) channel between the $m^{\mathrm{th}}$ vertex and the $k^{\mathrm{th}}$ user.
$T$	Completion time (in $s$) of the data collection and moving along the path.
$\beta$	Performance loss due to imperfect modulation and coding schemes.
$\eta$	Tag scattering efficiency.
$N_0$	Receiver noise power (in $\mathrm{Watt}$).
${\gamma }_k$	The communicaiton QoS target (in $\mathrm{bit}/\mathrm{Hz}$) at IoT user $k$.

Equations

$$\frac{1}{a}\sum _{m=1}^M\sum _{j=1}^M{W}_{m,j}{D}_{m,j}=\frac{\mathrm{Tr}({\mathbf{D}}^T\mathbf{W})}{a}.$$

$$E=\left(\frac{{\alpha }_1}{a}+{\alpha }_2\right)\mathrm{Tr}({\mathbf{D}}^T\mathbf{W}),$$

$$R_{k,m}={\log }_2\left(1+v_m\cdot \frac{\beta \eta |g_{k,m}{|}^2|h_{k,m}{|}^2{p}_{k,m}}{N_0}\right),$$

$$\sum _{m=1}^M{t}_{k,m}\cdot {\log }_2\left(1+v_m\cdot \frac{\beta \eta |g_{k,m}{|}^2|h_{k,m}{|}^2{p}_{k,m}}{N_0}\right)\ge {\gamma }_k,$$

$$\sum _{j=1}^M{W}_{m,j}=v_m,\sum _{j=1}^M{W}_{j,m}=v_m,\forall m=1,\cdots ,M.$$

$${\lambda }_m-{\lambda }_j+\left(\sum _{l=1}^M{v}_l-1\right)W_{m,j}+\left(\sum _{l=1}^M{v}_l-3\right)W_{j,m}$$

$$\le \sum _{l=1}^M{v}_l-2+J\left(2-v_m-v_j\right),\forall m,j\ge 2,m\ne j,$$

$$v_m\le {\lambda }_m\le \left(\sum _{l=1}^M{v}_l-1\right)v_m,\forall m\ge 2,$$

$$\mathrm{P}1:$$

$$+(2-\mu )\sum _{m=1}^M\sum _{k=1}^K{t}_{k,m}{p}_{k,m}$$

$$s.t.$$

$$\ge {\gamma }_k,\forall k,$$

$$\frac{1}{a}\mathrm{Tr}({\mathbf{D}}^T\mathbf{W})+\sum _{m=1}^M\sum _{k=1}^K{t}_{k,m}\le T,$$

$$\sum _{j=1}^M{W}_{m,j}=v_m,\sum _{j=1}^M{W}_{j,m}=v_m,\forall m,$$

$${\lambda }_m-{\lambda }_j+\left(\sum _{l=1}^M{v}_l-1\right)W_{m,j}+\left(\sum _{l=1}^M{v}_l-3\right)W_{j,m}$$

$$\le \sum _{l=1}^M{v}_l-2+J\left(2-v_m-v_j\right),$$

$$\forall m,j\ge 2,m\ne j,$$

$$v_m\le {\lambda }_m\le \left(\sum _{l=1}^M{v}_l-1\right)v_m,\forall m\ge 2,$$

$$W_{m,j}\in \{0,1\},\forall m,j,W_{m,m}=0,\forall m,$$

$$v_1=1,v_m\in \{0,1\},\forall m\ge 2,$$

$$(1-v_m)\cdot t_{k,m}=0,\forall k,m,$$

$$t_{k,m}\ge 0,p_{k,m}\ge 0,\forall k,m,$$

$$A_{k,m}=\frac{\beta \eta |g_{k,m}{|}^2|h_{k,m}{|}^2}{N_0},$$

$$t_{k,m}=\{\begin{array}{rlrl}& s_k,& \mathrm{if}m={\mathrm{argmax}}_{l\in \mathcal{V}}{v}_l{A}_{k,l}& \\ & 0,& \mathrm{if}m\ne {\mathrm{argmax}}_{l\in \mathcal{V}}{v}_l{A}_{k,l}& \end{array},$$

$$p_{k,m}=\{\begin{array}{rlrl}& q_k,& \mathrm{if}m={\mathrm{argmax}}_{l\in \mathcal{V}}{v}_l{A}_{k,l}& \\ & 0,& \mathrm{if}m\ne {\mathrm{argmax}}_{l\in \mathcal{V}}{v}_l{A}_{k,l}& \end{array},$$

$$\mathrm{P}2:{\min }_{\begin{array}{c}\mathbf{v},\mathbf{W},\{{\lambda }_m\}\\ \{s_k>0,q_k>0\}\end{array}}$$

$$+(2-\mu )\sum _{k=1}^K{s}_k{q}_k$$

$$s.t.$$

$$\ge {\gamma }_k,\forall k,$$

$$\frac{\mathrm{Tr}({\mathbf{D}}^T\mathbf{W})}{a}+\sum _{k=1}^K{s}_k\le T,$$

$$\text{eqref}P1c-\text{eqref}vertex.$$

$$q_k=\frac{1}{B_k(\tilde{\mathbf{v}})}\left(2^{{\gamma }_k/s_k}-1\right).$$

$$\mathrm{P}3:{\min }_{\begin{array}{c}\{s_k>0\}\end{array}}$$

$$+(2-\mu )\sum _{k=1}^K\frac{{\gamma }_k}{B_k(\tilde{\mathbf{v}})}\Theta \left(\frac{s_k}{{\gamma }_k}\right)$$

$$s.t.$$

$$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\eqref{P1c}-\eqref{edge}\Big{\}},$$

$${\max }_{\begin{array}{c}\mathbf{W},\{{\lambda }_m\}\end{array}}T-\frac{\mathrm{Tr}({\mathbf{D}}^T\mathbf{W})}{a},s.t.\text{eqref}P1c-\text{eqref}edge.$$

$${\rm Y}(\tilde{\mathbf{v}}):=T-\frac{\mathrm{Tr}({\mathbf{D}}^T\hat{\mathbf{W}})}{a}.$$

$$\mathcal{L}\left(\{s_k\},\rho \right)=$$

$$+(2-\mu )\sum _{k=1}^K\frac{{\gamma }_k}{B_k(\tilde{\mathbf{v}})}\Theta \left(\frac{s_k}{{\gamma }_k}\right)$$

$$+\rho \left(\sum _{k=1}^K{s}_k-{\rm Y}(\tilde{\mathbf{v}})\right).$$

$$-\nabla \Theta \left(\frac{{\hat{s}}_k}{{\gamma }_k}\right)=\frac{B_k(\tilde{\mathbf{v}})\hat{\rho }}{2-\mu },$$

$$-\nabla \Theta (x)=1+\ln 2\cdot 2^{1/x}/x-2^{1/x},$$

Full text

I Introduction

With a wide range of commercial and industrial applications, Internet of Things (IoT) market is continuously growing [1], 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 [2], 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 [3, 4, 5, 6, 7]. Unfortunately, due to the round-trip path-loss, the transmission range of backscatter communication is limited [8, 9, 10]. This can be seen from a recent prototype in [3], where the wirelessly powered backscatter communication only supports a range of $1$ meter at the data-rate of $1$ $\mathrm{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 [11]. 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 [12, 13, 14, 15]. 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 [16, 17, 18], 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 [19, 20, 21, 22, 23]. 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 $|\cdot|$ is the cardinality of a set. We use $(a_{1},a_{2},\cdots)$ to represent a sequence and $[a_{1},a_{2},\cdots]^{T}$ to represent a column vector, with $(\cdot)^{T}$ being the transpose operator. The operators $\textrm{Tr}(\cdot)$ and $(\cdot)^{-1}$ take the trace and the inverse of a matrix, respectively. Finally, $\mathbb{E}(\cdot)$ represents the expectation of a random variable.

II System Model

II-A 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 $(\mathcal{V},\mathcal{E})$ as shown in Fig. 1, where $\mathcal{V}=\{1,\cdots,M\}$ is the set of $M$ vertices representing the possible stopping points, and $\mathcal{E}$ is the set of directed edges representing the allowed movement paths [24]. To quantify the path length, a matrix $\mathbf{D}=[D_{1,1},\cdots,D_{1,M};\cdots;D_{M,1},\cdots,D_{M,M}]\in\mathbb{R}^{M\times 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}=+\infty$ [24]. Notice that $\mathbf{D}$ 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 $\mathcal{Q}=(y_{1},y_{2},\cdots,y_{Q})$ where $y_{j}\in\mathcal{V}$ for $j=1,\cdots,Q$ and $(y_{j},y_{j+1})\in\mathcal{E}$ for $j=1,\cdots,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 [16]. 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},\cdots,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\in\mathcal{V}$ [16, 25]. Thus this scenario can be represented by an extended graph with one more vertex and an extended $\mathbf{D}$ with dimension $(M+1)\times(M+1)$.

Correspondingly, we define the selection variable $\mathbf{v}=[v_{1},\cdots,v_{M}]^{T}\in\{0,1\}^{M}$, where $v_{m}=1$ if the vertex $m$ appears in the path $\mathcal{Q}$ and $v_{m}=0$ otherwise. Furthermore, we define a matrix $\mathbf{W}=[W_{1,1},\cdots,W_{1,M};\cdots;W_{M,1},\cdots,W_{M,M}]\in\{0,1\}^{M\times M}$, with $W_{y_{j},y_{j+1}}=1$ for all $j=1,\cdots,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 $\mathcal{Q}$ is

[TABLE]

Furthermore, since the total motion energy $E$ of the UGV is proportional to the total motion time [11, 12, 13], the motion energy can be expressed in the form of

[TABLE]

where $\alpha_{1}$ and $\alpha_{2}$ are parameters of the model (e.g., for a Pioneer 3DX robot in Fig. 1, $\alpha_{1}=0.29$ and $\alpha_{2}=7.4$ [11, Sec. IV-C]).

II-B Backscatter Data Collection Model

Based on the mobility model, the UGV moves along the selected path $\mathcal{Q}$ 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 communication111When user $k$ adapts the variable impedance for modulating the backscattered waveform with information bits, other users keep silent to avoid collision [6]. [26, 27]. 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}\neq 0$, the RF source at the UGV transmits a symbol $x_{k,m}\in\mathbb{C}$ with $\mathbb{E}[|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 $\eta|g_{k,m}|^{2}|h_{k,m}|^{2}p_{k,m}/N_{0}$, where $h_{k,m}\in\mathbb{C}$ is the downlink channel from the UGV to user $k$, $g_{k,m}\in\mathbb{C}$ 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 [26, 27]). Furthermore, $\eta$ is the tag scattering efficiency determined by the load impedance $Z_{L}$ and the antenna impedance $Z_{A}$ [28]. 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}\neq 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, $\eta=|(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

[TABLE]

where $\beta$ is the performance loss due to imperfect modulation and coding schemes in backscatter communication [29]. For example, in bistatic backscatter communication with frequency shift keying, $\beta=0.5$ [29]. On the other hand, in ambient backscatter communication with on-off keying, $\beta$ is obtained by fitting $\mathrm{log}_{2}\left(1+\beta x\right)$ to $1-\mathbb{Q}\left(\sqrt{x}\right)$ [30], where $\mathbb{Q}\left(x\right)=1/\sqrt{2\pi}\int_{x}^{\infty}\mathrm{exp}\left(-u^{2}/2\right)\mathrm{d}u$ 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 [31] 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) [32]. Then based on the probabilistic framework in [32], 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 $\eta$ to facilitate the analysis [8, 9, 29, 33]. On the other hand, if $\eta$ is not a constant, according to [34], it is possible to compute the range of $\eta$. Then, the worst case approach, which replaces $\eta$ in (3) with its lower bound $\eta_{\mathrm{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 $\eta$ as a nonlinear function of the incident power is an important future work.

III 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 $\mathbf{v}$ and $\mathbf{W}$) and designing the stopping time $\{t_{k,m}\}$ and transmit power $\{p_{k,m}\}$. In particular, the data collection QoS requirement of the $k^{\mathrm{th}}$ IoT device can be described by

[TABLE]

where $\gamma_{k}>0$ (in $\mathrm{bit/Hz}$) is the amount of data to be collected from user $k$.

Notice that the variables $\mathbf{v}$ and $\mathbf{W}$ are dependent since $v_{m}=0$ implies $W_{m,j}=W_{j,m}=0$ for any $j\in\mathcal{V}$. On the other hand, the UGV would visit the vertex with $v_{m}=1$, making $\sum_{j=1}^{M}W_{m,j}=\sum_{j=1}^{M}W_{j,m}=1$. Combining the above two cases, we have

[TABLE]

Furthermore, since the path must be connected, the following subtour elimination constraints are required to eliminate disjointed sub-tours [25]:

[TABLE]

where $\{\lambda_{m}\}$ are slack variables to guarantee a connected path, and $\sum_{l=1}^{M}v_{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 $\left(\alpha_{1}/a+\alpha_{2}\right)\mathrm{Tr}(\mathbf{D}^{T}\mathbf{W})$ and communication energy $\sum_{m=1}^{M}\sum_{k=1}^{K}t_{k,m}p_{k,m}$, the joint mobility management and power allocation problem of the data collection system is formulated as222The considered system adopts semi-passive backscatter communication, where the local circuits at IoT devices are powered by their own batteries [5]. 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 [35, 36, 37, 38, 39].:

[TABLE]

where (7b) is for constraining the operation (including moving and data collection) to be completed within $T$ seconds, and (7h) is for constraining the stopping time to be zero if the vertex is not visited. Notice that $0<\mu\leq 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 $\mu=1$. On the other hand, if we want to restrict the interference to other co-existing wireless systems, we might set $\mu<1$.

It can be seen from the constraint (7a) of $\mathrm{P}1$ 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 $\mathrm{P}1$ can concisely balance this energy trade-off.

Unfortunately, problem $\mathrm{P}1$ is nontrivial to solve due to the following reasons. Firstly, it is NP-hard, since it involves the integer constraints (7f)$-$(7g) [40]. 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 [40], where the reward of visiting each vertex is a constant.

IV Optimal Solution to $\mathrm{P}1$

Despite the optimization challenges, this section proposes an algorithm that theoretically obtains the optimal solution to $\mathrm{P}1$. The idea of this algorithm is to eliminate the variables $\mathbf{W},\{\lambda_{m}\},\{t_{k,m},p_{k,m}\}$ so as to transform $\mathrm{P}1$ into an equivalent problem only related to $\mathbf{v}$. 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 $\mathrm{P}1$ will be first discussed, which helps in reducing the dimension of $\mathrm{P}1$.

IV-A Optimality Condition

To address the challenges for solving $\mathrm{P}1$, we first establish the optimality condition of $\mathrm{P}1$. In particular, by defining

[TABLE]

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

Proposition 1**.**

The optimal $\{\mathbf{v}^{*},t^{*}_{k,m},p_{k,m}^{*}\}$ to $\mathrm{P}1$ satisfies:

(i) If $m=\mathop{\mathrm{argmax}}_{l\in\mathcal{V}}~{}v_{l}^{*}A_{k,l}$, then $t^{*}_{k,m}\neq 0$; otherwise $t^{*}_{k,m}=0$.

(ii) If $t_{k,m}^{*}\neq 0$, then $p^{*}_{k,m}\neq 0$.

Proposition 1 indicates that the UGV only needs to allocate time to user $k$ at a single vertex, which is given by $m=\mathrm{argmax}_{l\in\mathcal{V}}~{}v^{*}_{l}A_{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

[TABLE]

where $s_{k}>0$, without changing the optimal solution to $\mathrm{P}1$. Correspondingly, the transmit power $\{p_{k,m}\}$ can be set to333If $m\neq\mathop{\mathrm{argmax}}_{l\in\mathcal{V}}~{}v_{l}A_{k,l}$, we have $t_{k,m}=0$. Thus $t_{k,m}p_{k,m}=0$ in the objective of $\mathrm{P}1$ and $t_{k,m}\mathrm{log}_{2}\left(1+v_{m}A_{k,m}p_{k,m}\right)=0$ in (7a), meaning that $p_{k,m}$ would not participate in problem $\mathrm{P}1$. As a result, we can set $p_{k,m}=0$.

[TABLE]

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 (9) and (10) into $\mathrm{P}1$, problem $\mathrm{P}1$ is transformed into

[TABLE]

Notice that the constraint (7h) is dropped since (7h) is always satisfied when $\{t_{k,m}\}$ takes the form of (9).

The problem $\mathrm{P}2$ is still nontrivial to solve due to the nonlinear coupling between $\mathbf{v}$ and $\{\mathbf{W},s_{k},q_{k}\}$ as observed from the constraints (11a) and (11c). To resolve such coupling, a straightforward idea is to use alternating minimization for optimizing $\mathbf{v}$, $\mathbf{W}$ and $\{s_{k},q_{k}\}$ iteratively. However, due to the discrete nature of $\mathbf{v}$ and $\mathbf{W}$, 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 $\mathbf{W}$ and $\{s_{k},q_{k}\}$ to $\mathrm{P}2$ with fixed $\mathbf{v}$. By representing $\mathbf{W}$ and $\{s_{k},q_{k}\}$ as functions of $\mathbf{v}$, problem $\mathrm{P}2$ is simplified to an equivalent problem only involving $\mathbf{v}$. Then we will step further to find the optimal solution of vertex selection variable $\mathbf{v}$.

IV-B Optimal Solution of $\mathbf{W}$ and $\{s_{k},q_{k}\}$ with Fixed $\mathbf{v}$

When $\mathbf{v}=\widetilde{\mathbf{v}}$, where $\widetilde{\mathbf{v}}$ is any feasible solution to $\mathrm{P}2$, the constraint (7g) can be dropped since it only involves $\mathbf{v}$. On the other hand, the term $\mathop{\mathrm{max}}_{l\in\mathcal{V}}~{}\widetilde{v}_{l}A_{k,l}$ in constraint (11a) becomes a constant, and we denote it as $B_{k}(\widetilde{\mathbf{v}}):=\mathop{\mathrm{max}}_{l}~{}\widetilde{v}_{l}A_{k,l}$. Furthermore, it can be seen from the objective function of $\mathrm{P}2$ that $q_{k}$ is a variable to be minimized and $s_{k}\mathrm{log}_{2}\left[1+B_{k}(\widetilde{\mathbf{v}})q_{k}\right]$ is a strictly increasing function of $q_{k}$. As a result, the optimal solution of $q_{k}$ must activate the constraint (11a) of $\mathrm{P}2$, which leads to

[TABLE]

Putting (12) into $\mathrm{P}2$, it is proved in Appendix B that the optimal $\mathbf{W}^{*}$ and $\{s_{k}^{*}\}$ to $\mathrm{P}2$ must activate the constraint (11b), i.e., $\mathrm{Tr}(\mathbf{D}^{T}\mathbf{W}^{*})/a+\sum_{k=1}^{K}s_{k}^{*}=T$. Using this result, the quantity $\mathrm{Tr}(\mathbf{D}^{T}\mathbf{W})$ in the objective function of $\mathrm{P}2$ can be replaced by $a(T-\sum_{k=1}^{K}s_{k})$ without changing the problem, and the objective function would be independent of $\mathbf{W}$. Then problem $\mathrm{P}2$ is equivalently transformed into the following two-stage optimization problem (detailed procedure given in Appendix B):

[TABLE]

where $\Theta(x):=x(2^{1/x}-1)$.

To solve $\mathrm{P}3$, we first need to compute the right hand side of the constraint, which leads to the following problem:

[TABLE]

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

Denoting the optimal solution to the problem (13) as $\{\widehat{\mathbf{W}},\widehat{\lambda}_{m}\}$, the optimal objective value of the travelling salesman problem is given by

[TABLE]

Now, by putting the obtained $\Upsilon(\widetilde{\mathbf{v}})$ into $\mathrm{P}3$, the constraint of $\mathrm{P}3$ is written as $\sum_{k=1}^{K}s_{k}=\Upsilon(\widetilde{\mathbf{v}})$. Assigning a Lagrange multiplier $\rho$ to this constraint, the Lagrangian of $\mathrm{P}3$ is

[TABLE]

According to the first-order KKT condition $\frac{\partial\mathcal{L}}{\partial s_{k}}|_{s_{k}=\widehat{s}_{k}}=0$, the optimal $\{\widehat{s}_{k}\}$ and $\widehat{\rho}$ should together satisfy [41]:

[TABLE]

where $-\nabla\Theta(x)$ is the gradient of $-\Theta(x)$ and is given by

[TABLE]

with $x>0$. Moreover, as $-\nabla^{2}\Theta(x)=-\mathrm{ln}^{2}2\cdot 2^{1/x}/x^{3}<0$ and $-\nabla^{3}\Theta(x)=\mathrm{ln}^{3}2\cdot 2^{1/x}/x^{5}+3\mathrm{ln}^{2}2\cdot 2^{1/x}/x^{4}>0$, it can be seen that $-\nabla\Theta(x)$ is a strictly decreasing and convex function of $x$. As a result, there must exist a strictly decreasing and convex function $\Lambda(x)$ such that $\Lambda(-\nabla\Theta(x))=x$. That is, the function $\Lambda(x)$ is the inverse function of $-\nabla\Theta(x)$, and it can be numerically computed and stored as a look-up table, with its shape shown in Fig. 3. Applying $\Lambda(x)$ to both sides of (15), we have

[TABLE]

Notice that in (17), the only unknown is $\widehat{\rho}$, which should satisfy the equality constraint of $\mathrm{P}3$:

[TABLE]

Since $\Lambda(x)$ is a decreasing function, bisection search method can be used to find $\widehat{\rho}$ efficiently. In order to determine the bisection interval, the following proposition (proved in Appendix C) can be established.

Proposition 2**.**

The quantity $\widehat{\rho}$ in (18) is bounded as

[TABLE]

Once $\widehat{\rho}$ is obtained, we can put $\widehat{\rho}$ into (17) to get $\{\widehat{s}_{k}\}$. Further putting $\{\widehat{s}_{k}\}$ into (12), we have \widehat{q}_{k}=\left(2^{\gamma_{k}/\widehat{s}_{k}}-1\right)\Big{/}B_{k}(\widetilde{\mathbf{v}}).

IV-C Optimal Solution of $\mathbf{v}$

With path $\widehat{\mathbf{W}}$, transmit times $\{\widehat{s}_{k}\}$, and transmit powers $\{\widehat{q}_{k}\}$ derived in Section IV-B, the optimal objective value of $\mathrm{P}3$ (equivalently $\mathrm{P}2$ with $\mathbf{v}=\widetilde{\mathbf{v}}$) is given by

[TABLE]

Therefore, problem $\mathrm{P}2$ is re-written as

[TABLE]

To solve $\mathrm{P}4$, a naive way is to apply exhaustive search for $\mathbf{v}$. 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 $\mathbf{v}$, leading to significant reduction of the computational complexity compared to exhaustive search while guaranteeing the global optimality [42, 43].

In particular, we define the living pool as a set $\mathcal{Y}$ 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 $\mathbf{v}$ for $\mathrm{P}4$ is $\Omega=\left\{\mathbf{v}\in\{0,1\}^{M}:v_{1}=1\right\}$, the initial living pool can be set to $\mathcal{Y}=\{\mathcal{F}_{(1,0)},\mathcal{F}_{(1,1)}\}$ (notice that $\mathcal{Y}$ is a family of sets over $\Omega$), where $\mathcal{F}_{(1,0)}=\left\{\mathbf{v}\in\Omega:v_{2}=0\right\}$ and $\mathcal{F}_{(1,1)}=\left\{\mathbf{v}\in\Omega:v_{2}=1\right\}$ as shown in Fig. 4a. It can be seen that $\Omega=\mathcal{F}_{(1,0)}\bigcup\mathcal{F}_{(1,1)}$ and $\mathcal{F}_{(1,0)},\mathcal{F}_{(1,1)}\subset\Omega$.

At the beginning, the B$\&$B method computes a lower bound $\Psi(\mathcal{F}_{(1,0)})$ such that $\Psi(\mathcal{F}_{(1,0)})\leq\Xi(\widetilde{\mathbf{v}})$ for any $\widetilde{\mathbf{v}}\in\mathcal{F}_{(1,0)}$, where $\Xi(\widetilde{\mathbf{v}})$ is defined in (IV-C). Based on the bounding function value $\Psi(\mathcal{F}_{(1,0)})$, we consider the following three cases.

• (i)

If $\Psi(\mathcal{F}_{(1,0)})>I$, then the subset $\mathcal{F}_{(1,0)}$ can be discarded, since no feasible solution inside $\mathcal{F}_{(1,0)}$ leads to better objective than the incumbent. Therefore, we update $\mathcal{Y}\leftarrow\mathcal{Y}\setminus\mathcal{F}_{(1,0)}$. This case is shown in Fig. 4b.

• (ii)

If $\Psi(\mathcal{F}_{(1,0)})\leq I$ and $|\mathcal{F}_{(1,0)}|>1$, the possibility of a better solution in $\mathcal{F}_{(1,0)}$ cannot be ruled out. As a result, we need to branch on $\mathcal{F}_{(1,0)}$ and generates two subsets $\mathcal{F}_{(1,0,0)}=\left\{\mathbf{v}\in\mathcal{F}_{(1,0)}:v_{3}=0\right\}$ and $\mathcal{F}_{(1,0,1)}=\left\{\mathbf{v}\in\mathcal{F}_{(1,0)}:v_{3}=1\right\}$. By treating $\mathcal{F}_{(1,0,0)}$ and $\mathcal{F}_{(1,0,1)}$ as new subset nodes, we update $\mathcal{Y}\leftarrow\mathcal{Y}\setminus\mathcal{F}_{(1,0)}\cup\{\mathcal{F}_{(1,0,0)},\mathcal{F}_{(1,0,1)}\}$. This case is shown in Fig. 4c.

• (iii)

Otherwise, we must have $|\mathcal{F}_{(1,0)}|=1$. Denoting the unique element in $|\mathcal{F}_{(1,0)}|$ as $\widetilde{\mathbf{v}}$, we compute $\Xi(\widetilde{\mathbf{v}})$. If $\Xi(\widetilde{\mathbf{v}})>I$, $\mathcal{F}_{(1,0)}$ is discarded. On the other hand, if $\Xi(\widetilde{\mathbf{v}})\leq I$, the incumbent is updated as $I\leftarrow\Xi(\widetilde{\mathbf{v}})$ and the current best solution is updated as $\mathbf{v}^{\diamond}\leftarrow\widetilde{\mathbf{v}}$. In both cases, since there is no more element in $\mathcal{F}_{(1,0)}$ to be evaluated, we can update $\mathcal{Y}\leftarrow\mathcal{Y}\setminus\mathcal{F}_{(1,0)}$.

After $\mathcal{F}_{(1,0)}$ is evaluated, we then evaluate $\mathcal{F}_{(1,1)},\mathcal{F}_{(1,0,0)},\mathcal{F}_{(1,0,1)},\cdots$ by repeating the above procedure until the living pool $\mathcal{Y}$ becomes empty.

For the above B$\&$B method, the key step is to derive the bounding function $\Psi$. In particular, to guarantee that the proposed B$\&$B method finds the optimal solution of $\mathbf{v}$ to $\mathrm{P}4$, the bounding function needs to satisfy [42]:

[TABLE]

where $(z_{1},\cdots,z_{N})$ are the values assigned to $(v_{1},\cdots,v_{N})$ and $N$ is the number of fixed elements in $\mathbf{v}$. To this end, consider the following bounding function $\Psi$:

[TABLE]

where $\mathbf{z}=[z_{1},\cdots,z_{N}]^{T}\in\{0,1\}^{N}$ (corresponding to the sequence $(z_{1},\cdots,z_{N})$), and $\delta$ is the solution to

[TABLE]

The function

[TABLE]

represents a bipartite matching problem, which can be numerically computed via the Hungarian algorithm [44]. It is proved in Appendix D that (23) satisfies the property (22). As a result, by applying the bounding function $\Psi$ in (23), the proposed B$\&$B method is guaranteed to obtain the optimal solution of $\mathbf{v}$ to problem $\mathrm{P}4$ (equivalently $\mathrm{P}2$).

IV-D Summary of Algorithm and Complexity Analysis

Since the B$\&$B algorithm finds the optimal solution of $\mathbf{v}$ to $\mathrm{P}2$, and the optimal solution of $\mathbf{W}$ and $\{s_{k},q_{k}\}$ with fixed $\mathbf{v}$ can be computed according to Section IV-B, the entire algorithm for computing the optimal solution to $\mathrm{P}2$ (equivalently $\mathrm{P}1$) is summarized in Algorithm 1.

In terms of computational complexity, computing $\Phi\left(\mathbf{z}\right)$ in (25) via the Hungarian algorithm requires a complexity of $O\left((2M)^{3}\right)$ [44]. On the other hand, computing $\delta$ would require bisection search in solving (24) and the number of iterations is given by $\mathrm{log}_{2}\left(\frac{C}{\epsilon}\right)$ [45], where $C$ is the length of the initial searching interval given by (19) and $\epsilon$ is the target accuracy. Therefore, computing the bounding function $\Psi$ requires a complexity of $O\left((2M)^{3}+K\mathrm{log}_{2}\left(\frac{C}{\epsilon}\right)\right)$. Finally, computing $\Xi\left(\widetilde{\mathbf{v}}\right)$ for a fixed $\mathbf{v}=\widetilde{\mathbf{v}}$ would involve the travelling salesman problem, which requires a complexity of $O\left((M-1)^{2}\cdot 2^{M-1}\right)$ in the worst case [46].

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

[TABLE]

where $X_{1}$ is the number of times for computing $\Psi$ and $X_{2}$ is the number of times for computing $\Xi\left(\widetilde{\mathbf{v}}\right)$. In contrast, the computational complexity for exhaustive search of $\mathbf{v}$ is $2^{M-1}\cdot(M-1)^{2}\cdot 2^{M-1}$. As $(2M)^{3}+K\mathrm{log}_{2}(C/\epsilon)$ is much smaller than $(M-1)^{2}\cdot 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 $\mathrm{P}1$. However, since we are minimizing the objective function of $\mathrm{P}1$, a smaller initial incumbent $I$ could help in reducing the size of the living pool [47], 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}\leq P_{\mathrm{max}}$ (with $P_{\mathrm{max}}$ being the upper limit of transmit power) for all $k$ to problem $\mathrm{P}3$. In such a case, the proposed Algorithm 1 is still applicable if we modify the objective function $\Xi(\mathbf{v})$ in problem $\mathrm{P}4$ into

[TABLE]

and the bounding function $\Psi$ in (23) into

[TABLE]

V 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 [48, 49]. More specifically, we start from a feasible solution of $\mathbf{v}$ (e.g., $\mathbf{v}^{[0]}=[1,0,\cdots,0]^{T}$), and randomly selects a candidate solution $\mathbf{v}^{\prime}$ from the neighborhood $\mathcal{N}(\mathbf{v}^{[0]})$. Since a natural neighborhood operator for binary optimization is to flip the value of $\{v_{m}\}$, $\mathcal{N}(\mathbf{v}^{[0]})$ can be set to

[TABLE]

where $L\geq 1$ is the size of neighborhood. It can be seen that $\mathcal{N}(\mathbf{v}^{[0]})$ is a subset of the entire feasible space $\Omega$ and containing solutions “close” to $\mathbf{v}^{[0]}$.

With the neighborhood $\mathcal{N}(\mathbf{v}^{[0]})$ defined above and the choice of $\mathbf{v}$ fixed to $\mathbf{v}=\mathbf{v}^{\prime}$, we consider two cases.

• (i)

If $\Xi(\mathbf{v}^{\prime})\leq\Xi(\mathbf{v}^{[0]})$, we update $\mathbf{v}^{[1]}\leftarrow\mathbf{v}^{\prime}$. By treating $\mathbf{v}^{[1]}$ as a new feasible solution, we can construct the next neighborhood $\mathcal{N}(\mathbf{v}^{[1]})$.

• (ii)

If $\Xi(\mathbf{v}^{\prime})>\Xi(\mathbf{v}^{[0]})$, we find another point within the neighborhood $\mathcal{N}(\mathbf{v}^{[0]})$ until $\Xi(\mathbf{v}^{\prime})\leq\Xi(\mathbf{v}^{[0]})$.

The above procedure is repeated to generate a sequence of $\{\mathbf{v}^{[1]},\mathbf{v}^{[2]},\cdots\}$ and the converged point is guaranteed to be a local optimal solution to $\mathrm{P}1$ [40]. 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 $\overline{\mathrm{Iter}}$. Denoting the solution after $\overline{\mathrm{Iter}}$ iterations as $\mathbf{v}^{\diamond}$, we can set the initial incumbent as $I=\Xi(\mathbf{v}^{\diamond})$. The entire procedure to generate an initial incumbent is summarized in Algorithm 2, and the complexity of executing Algorithm 2 is $\overline{\mathrm{Iter}}\cdot(M-1)^{2}\cdot 2^{M-1}$.

VI 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 $\eta=0.78$ (corresponding to $1.1~{}\mathrm{dB}$ loss [29]), the performance loss due to imperfect modulation is $\beta=0.5$ [29], and the weighting factor $\mu=1$. Within the time budget $T=500~{}\mathrm{s}$, the data collection targets $\gamma_{k}\sim\mathcal{U}(1,2)$ in the unit of $\mathrm{bit/Hz}$ are requested by $K=50$ IoT users (corresponding to a spectral efficiency of $K\gamma_{k}/T=0.1\sim 0.2~{}\mathrm{bps/Hz}$ [1]), where $\mathcal{U}(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~{}\mathrm{m}\times 20~{}\mathrm{m}=400~{}\mathrm{m}^{2}$ square area, which is a typical size for smart warehouses. Inside this map, $K=50$ IoT users444Notice 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 $\varrho_{k,m}=\varrho_{0}\cdot(\frac{d_{k,m}}{d_{0}})^{-2.5}$ is adopted [50], where $d_{k,m}$ is the distance from user $k$ to the stopping point $m$, and $\varrho_{0}=10^{-3}$ is the path-loss at $d_{0}=1~{}\mathrm{m}$. Based on the path-loss model, channels $g_{k,m}$ and $h_{k,m}$ are generated according to $\mathcal{CN}(0,\varrho_{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 $\mathcal{Y}$ versus the number of iterations (represented by $\mathrm{Iter}$ in Algorithm 1) when the receiver noise power $N_{0}=-95~{}\mathrm{dBm}$ (corresponding to power spectral density $-145~{}\mathrm{dBm/Hz}$ with $100~{}\mathrm{kHz}$ bandwidth [1]). It can be seen that even with a naive initial incumbent obtained by setting $\mathbf{v}=[1,0,\cdots,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 $\mathrm{P}1$ with $\mathbf{v}=[1,0,\cdots,0]^{T}$) and the scheme with full path visiting all vertices (i.e., the optimal solution to $\mathrm{P}1$ with $\mathbf{v}=\mathbf{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~{}\mathrm{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~{}\mathrm{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$ $\mathrm{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$ $\mathrm{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~{}\mathrm{W}$ (occurred at noise power $-60~{}\mathrm{dBm}$). Translating this number to the received power density at $1~{}\mathrm{m}$ gives $2.45~{}\mathrm{W/m^{2}}$, which is within the requirement ($<10~{}\mathrm{W/m^{2}}$) set by the IEEE standard C95.1-2005 [2, Remark 4].

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 $\eta$ could also impact the path. To see this, the case of $K=50$ with $M=12$ at noise power $N_{0}=-80~{}\mathrm{dBm}$ is simulated, and the paths for $(T,\eta)=(500~{}\mathrm{s},0.78)$, $(T,\eta)=(20~{}\mathrm{s},0.78)$, and $(T,\eta)=(500~{}\mathrm{s},0.1)$ are compared in Fig. 7. It can be seen from Fig. 7 that the path for $(T,\eta)=(20~{}\mathrm{s},0.78)$ involves a smaller moving distance than that for $(T,\eta)=(500~{}\mathrm{s},0.78)$. This is because a smaller $T$ would restrict the constraint (7b) of $\mathrm{P}1$, which forces the UGV to reduce its motion time and moving distance. On the other hand, the path for $(T,\eta)=(500~{}\mathrm{s},0.1)$ involves a larger moving distance than that for $(T,\eta)=(500~{}\mathrm{s},0.78)$, since a smaller $\eta$ 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 $\gamma_{1}=\gamma_{2}=\cdots=\gamma_{K}$ in $\mathrm{bit/Hz}$ at noise power $N_{0}=-90$ $\mathrm{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 $\mu$ varies, the case of $K=50$ with $M=12$ at $N_{0}=-90~{}\mathrm{dBm}$ is simulated and the result is shown in Fig. 9. It can be seen that if $\mu=0$, the proposed Algorithm 1 has the same performance as that of the full path. This is because $\mu=0$ would make the motion energy disappear in the objective function of $\mathrm{P}1$, and the best strategy is to allow the UGV to visit all vertices. However, even with a slight increase in $\mu$, the proposed Algorithm 1 would outperform other benchmark schemes. Finally, it can be seen that if $\mu$ increases, the motion energy decreases but the communication energy increases. Therefore, $\mu$ 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 $\mu=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 [22, 23].

Finally, we analyze the impact of the number of vertices $M$ on the energy consumption. In particular, the case of $K=50$ with $M\in\{3,6,9,12\}$ at noise power $N_{0}=-90~{}\mathrm{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 [51].

VII 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.

Appendix A Proof of Proposition 1

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

A-A Proof of Part (ii)

To prove part (ii) by contradiction, consider an optimal solution $\{t_{k,m}^{*},p_{k,m}^{*},\mathbf{v}^{*},\mathbf{W}^{*},\lambda_{m}^{*}\}$ to $\mathrm{P}1$ with a particular $(i,j)$ such that $t_{i,j}^{*}=\widetilde{t}\neq 0$. Assume that the corresponding $p_{i,j}^{*}=0$.

Since $\{t_{k,m}^{*},p_{k,m}^{*},\mathbf{v}^{*},\mathbf{W}^{*},\lambda_{m}^{*}\}$ is optimal to $\mathrm{P}1$, it must satisfy (7a) of $\mathrm{P}1$, i.e., there must exist some $n\neq j$ such that $t_{i,n}^{*}=t^{\prime}\neq 0$ and $p_{i,n}^{*}=p^{\prime}\neq 0$. Now, consider a related problem of $\mathrm{P}1$ by fixing all the variables to their optimal values except for $(t_{i,j},t_{i,n},p_{i,j},p_{i,n})$:

[TABLE]

where

[TABLE]

As $\{t_{k,m}^{*},p_{k,m}^{*},\mathbf{v}^{*},\mathbf{W}^{*}\}$ is optimal to $\mathrm{P}1$, it can be seen that

[TABLE]

Therefore, $(\widetilde{t},t^{\prime},0,p^{\prime})$ should satisfy the constraints of (30), which leads to

[TABLE]

Furthermore, by Jensen’s inequality, we have

[TABLE]

where the strict inequality is due to $\widetilde{t},t^{\prime},p^{\prime}\neq 0$. Adding $C_{1}$ to both sides of (33), and combining the first inequality of (32), we have $C_{1}+(\widetilde{t}+t^{\prime})\mathrm{log}_{2}\left(1+v_{n}^{*}A_{i,n}\cdot\frac{t^{\prime}}{\widetilde{t}+t^{\prime}}p^{\prime}\right)>\gamma_{i}$. Comparing this result and the second inequality of (32) to the constraints of (30), it can be seen that $(t_{i,j},t_{i,n},p_{i,j},p_{i,n})=(0,\widetilde{t}+t^{\prime},0,\frac{t^{\prime}}{\widetilde{t}+t^{\prime}}p^{\prime}-\Delta p)$ is feasible for (30) under sufficiently small $\Delta p>0$. Putting $(0,\widetilde{t}+t^{\prime},0,\frac{t^{\prime}}{\widetilde{t}+t^{\prime}}p^{\prime}-\Delta p)$ and $(\widetilde{t},t^{\prime},0,p^{\prime})$ from (31) into the objective function of (30), we obtain $t^{\prime}p^{\prime}-(\widetilde{t}+t^{\prime})\Delta p$ and $t^{\prime}p^{\prime}$, respectively. Further due to $t^{\prime}p^{\prime}-(\widetilde{t}+t^{\prime})\Delta p<t^{\prime}p^{\prime}$, it is clear that $(\widetilde{t},t^{\prime},0,p^{\prime})$ cannot be optimal to (30). This contradicts to (31). Therefore, $p_{i,j}^{*}\neq 0$.

A-B 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}^{*},\mathbf{v}^{*},\mathbf{W}^{*},\lambda_{m}^{*}\}$ to $\mathrm{P}1$. Assume that there exists some user $i$ such that $t_{i,j}^{*}=0$ at vertex $j=\mathop{\mathrm{argmax}}_{l\in\mathcal{V}}~{}v_{l}^{*}A_{i,l}$.

Since $\{t_{k,m}^{*},p_{k,m}^{*},\mathbf{v}^{*},\mathbf{W}^{*},\lambda_{k}^{*}\}$ is optimal to $\mathrm{P}1$, it must satisfy (7a) of $\mathrm{P}1$, i.e., there must exist some $n\neq j$ such that $t_{i,n}^{*}=t^{\prime}\neq 0$ and $p_{i,n}^{*}=p^{\prime}\neq 0$. Furthermore, under $\{t_{k,m}^{*},p_{k,m}^{*},\mathbf{v}^{*},\mathbf{W}^{*},\lambda_{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^{\prime},p_{i,j}^{*},p^{\prime})$ is optimal to (30). As a result, $(0,t^{\prime},p_{i,j}^{*},p^{\prime})$ must satisfy the constraints of (30), i.e.,

[TABLE]

which can be re-written as

[TABLE]

Comparing (35) with the constraints of (30), it can be seen that $(t_{i,j},t_{i,n},p_{i,j},p_{i,n})=(t^{\prime},0,\frac{v_{n}^{*}A_{i,n}}{v_{j}^{*}A_{i,j}}p^{\prime},0)$ is also feasible for (30). Putting $(0,t^{\prime},p_{i,j}^{*},p^{\prime})$ and $(t^{\prime},0,\frac{v_{n}^{*}A_{i,n}}{v_{j}^{*}A_{i,j}}p^{\prime},0)$ into the objective function of (30), we obtain $p^{\prime}t^{\prime}$ and $\frac{v_{n}^{*}A_{i,n}}{v_{j}^{*}A_{i,j}}p^{\prime}t^{\prime}$, respectively. Further due to $(0,t^{\prime},p_{i,j}^{*},p^{\prime})$ being optimal to (30), we have $p^{\prime}t^{\prime}\leq\frac{v_{n}^{*}A_{i,n}}{v_{j}^{*}A_{i,j}}p^{\prime}t^{\prime}$, which leads to $v_{j}^{*}A_{i,j}\leq v_{n}^{*}A_{i,n}$. As $A_{i,j}$ and $A_{i,n}$ are not equal almost surely, and $v_{j}^{*},v_{n}^{*}\in\{0,1\}$, the equality sign in $v_{j}^{*}A_{i,j}\leq 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=\mathop{\mathrm{argmax}}_{l\in\mathcal{V}}~{}v_{l}^{*}A_{i,l}$. Therefore, $t_{i,j}^{*}\neq 0$ at vertex $j=\mathop{\mathrm{argmax}}_{l\in\mathcal{V}}~{}v_{l}^{*}A_{i,l}$.

A-C 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}\neq 0$ at vertex $n\neq\mathop{\mathrm{argmax}}_{l\in\mathcal{V}}~{}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}\neq 0$ at vertex $j=\mathop{\mathrm{argmax}}_{l\in\mathcal{V}}~{}v_{l}^{*}A_{i,l}$. Correspondingly, by part (ii) of Proposition 1, we have $p^{*}_{i,n},p^{*}_{i,j}\neq 0$.

Now, the partial Lagrangian of $\mathrm{P}1$ with respect to $\{t_{k,m},p_{k,m}\}$ under fixed $\mathbf{v}=\mathbf{v}^{*}$ is

[TABLE]

where $\{\zeta_{k},\varphi,\chi_{k,m},\theta_{k,m},\xi_{k,m}\}$ are Lagrange multipliers. Since $\mathrm{P}1$ is convex in $\{t_{k,m}\}$ with $\{p_{k,m}\}$ fixed (vice versa), according to the Karush-Kuhn-Tucker condition [41], the optimal $\{t^{*}_{k,m},p^{*}_{k,m}\}$ and $\{\zeta_{k}^{*},\varphi^{*},\chi_{k,m}^{*},\theta_{k,m}^{*},\xi_{k,m}^{*}\}$ must satisfy

[TABLE]

Putting $t^{*}_{i,j},t^{*}_{i,n},p_{i,j}^{*},p_{i,n}^{*}\neq 0$ into (36a), we have

[TABLE]

Further putting $\xi^{*}_{i,n}=0$ from (37) into (36b), the following equation is obtained

[TABLE]

Substituting (38), $\theta^{*}_{i,n}=0$ from (37), and $v_{n}^{*}=1$ (due to (7h) and $t^{*}_{i,n}\neq 0$) into (36c), equation (36c) is reformulated as $F(v_{n}^{*}A_{i,n})=T\varphi^{*}$, where

[TABLE]

with $x\neq 0$. Notice that due to $v_{n}^{*}A_{i,n}\neq 0$ (as $v_{n}^{*}=1$), $F(v_{n}^{*}A_{i,n})$ is well-defined. Similarly, by using $\xi^{*}_{i,j}=\theta^{*}_{i,j}=0$ from (37), we obtain $F(v_{j}^{*}A_{i,j})=T\varphi^{*}$. Therefore,

[TABLE]

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

[TABLE]

where the second equality is obtained from (38). Since $p^{*}_{i,n}>0$, it is clear that $\nabla_{x}F(x)>0$ for any $x\in[v^{*}_{n}A_{i,n},v^{*}_{j}A_{i,j}]$ and $F(x)$ is a strictly increasing function of $x$ over this interval. Combining the result from (40), we have $v^{*}_{n}A_{i,n}=v_{j}^{*}A_{i,j}=\mathop{\mathrm{max}}_{l\in\mathcal{V}}~{}v^{*}_{l}A_{i,l}$. This contradicts to $n\neq\mathop{\mathrm{argmax}}_{l\in\mathcal{V}}~{}v_{l}^{*}A_{i,l}$. Therefore, $t^{*}_{i,n}=0$ at vertex $n\neq\mathop{\mathrm{argmax}}_{l\in\mathcal{V}}~{}v_{l}^{*}A_{i,l}$.

Appendix B Transformation from $\mathrm{P}2$ to $\mathrm{P}3$ with Fixed $\mathbf{v}=\widetilde{\mathbf{v}}$

With $q_{k}$ in (12), the first constraint of $\mathrm{P}2$ is always satisfied. Therefore, we can re-write $\mathrm{P}2$ as

[TABLE]

Since $s_{k}$ is not involved in $\eqref{P1c}-\eqref{edge}$, the partial Lagrangian of the problem (42) with respect to $\{s_{k}\}$ is given by

[TABLE]

where $\varphi\geq 0$ and $\theta_{k}\geq 0$ are Lagrange multipliers. Since problem (42) is convex in $\{s_{k}\}$, we must have $\frac{\partial\mathcal{L}}{\partial s_{k}}=0$ according to the KKT condition, and the optimal $\{\widehat{s}_{k}\}$ and $(\widehat{\varphi},\widehat{\theta}_{k})$ should together satisfy

[TABLE]

Since the derivative of the left hand side of (43) with respect to $\widehat{s}_{k}$ is $-\frac{(2-\mu)}{B_{k}(\widetilde{\mathbf{v}})}\cdot\mathrm{ln}^{2}2\cdot 2^{\gamma_{k}/\widehat{s}_{k}}\cdot\gamma_{k}^{2}/(\widehat{s}_{k})^{3}<0,$ the left hand side function in (43) is a strictly decreasing function of $\widehat{s}_{k}$. Furthermore, due to

[TABLE]

it is clear that $\widehat{\varphi}-\widehat{\theta}_{k}>0$. As $\widehat{\theta}_{k}\widehat{s}_{k}=0$ from the complementary slackness condition, and due to $\widehat{s}_{k}>0$, we must have $\widehat{\theta}_{k}=0$, and therefore $\widehat{\varphi}>0$.

Finally, from $\widehat{\varphi}>0$ and the complementary slackness condition $\widehat{\varphi}\left[\frac{1}{a}\mathrm{Tr}(\mathbf{D}^{T}\widehat{\mathbf{W}})+\sum_{k=1}^{K}\widehat{s}_{k}-T\right]=0$, the equation $\mathrm{Tr}(\mathbf{D}^{T}\widehat{\mathbf{W}})=a\left(T-\sum_{k=1}^{K}\widehat{s}_{k}\right)$ holds. Substituting this result into $\mathrm{P}2$, $\mathrm{P}2$ becomes

[TABLE]

Now the objective function of (B) is independent of $\mathbf{W}$ 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 (B) is obtained when $s_{k}$ is maximized. Since $\sum_{k=1}^{K}s_{k}=T-\mathrm{Tr}(\mathbf{D}^{T}\mathbf{W})/a$, maximizing $T-\mathrm{Tr}(\mathbf{D}^{T}\mathbf{W})/a$ helps in enlarging the maximum values of $\{s_{k}\}$. Therefore, (B) is equivalently transformed into $\mathrm{P}3$.

Appendix C Proof of Proposition 2

We first prove the left part of equation (19). In particular, since $\Lambda(x)$ is a convex function (see Fig. 3), we must have

[TABLE]

where “$\leq$” is due to Jensen’s inequlity and “$=$” is due to (18). By further applying function $-\nabla\Theta$ to both sides of (46), and since $-\nabla\Theta$ is a strictly decreasing function, (46) becomes

[TABLE]

which immediately leads to the left part of (19).

Next, we prove the right part of equation (19). More specifically, since $B_{k}(\widetilde{\mathbf{v}})\geq\mathrm{min}_{l}~{}B_{l}(\widetilde{\mathbf{v}})$ and $\widehat{\rho}\geq 0$, we have

[TABLE]

Applying $\Lambda(x)$ to both sides of (47), and since $\Lambda(x)$ is a strictly decreasing function, (47) becomes

[TABLE]

Based on the above result and equation (18), it is clear that

[TABLE]

which is equivalent to

[TABLE]

This leads to the right part of (19), and the proof is completed.

Appendix D Proof of Lower Bound Property of $\Psi$ in (23)

To begin with, it is noticed that $\Xi\left(\widetilde{\mathbf{v}}\right)$ in (IV-C) is the optimal value of problem $\mathrm{P}3$. Therefore, to prove $\Psi(\mathcal{F}_{(z_{1},\cdots,z_{N})})\leq\Xi(\widetilde{\mathbf{v}})$ for any $\widetilde{\mathbf{v}}\in\mathcal{F}_{(z_{1},\cdots,z_{N})}$, we consider two relaxations for $\mathrm{P}3$ with $\widetilde{\mathbf{v}}\in\mathcal{F}_{(z_{1},\cdots,z_{N})}$.

(i) Based on the definition $B_{k}(\widetilde{\mathbf{v}}):=\mathop{\mathrm{max}}_{l}~{}\widetilde{v}_{l}A_{k,l}$, it can be seen that $B_{k}(\mathbf{v}_{1})\geq B_{k}(\mathbf{v}_{2})$ if $\mathbf{v}_{1}\succeq\mathbf{v}_{2}$, where $\succeq$ means “Pareto dominance”. Furthermore, since $[\mathbf{z}^{T},\mathbf{1}_{M-N}^{T}]^{T}\succeq\widetilde{\mathbf{v}}$ for any $\widetilde{\mathbf{v}}\in\mathcal{F}_{(z_{1},\cdots,z_{N})}$, we have $B_{k}\left([\mathbf{z}^{T},\mathbf{1}_{M-N}^{T}]^{T}\right)\geq B_{k}(\widetilde{\mathbf{v}})$ for any $\widetilde{\mathbf{v}}\in\mathcal{F}_{(z_{1},\cdots,z_{N})}$. Using this result, the objective function of $\mathrm{P}3$ can be lower bounded by $\mu\left(\alpha_{1}+\alpha_{2}a\right)\left(T-\sum_{k=1}^{K}s_{k}\right)+(2-\mu)\sum_{k=1}^{K}\frac{\gamma_{k}}{B_{k}\left([\mathbf{z}^{T},\mathbf{1}_{M-N}^{T}]^{T}\right)}\Theta\left(\frac{s_{k}}{\gamma_{k}}\right)$.

(ii) Notice that

[TABLE]

where the inequality is obtained by dropping constraints (7d) and (7e). Using (51), the constraint of $\mathrm{P}3$ can be relaxed into $\sum_{k=1}^{K}s_{k}\leq\Phi(\mathbf{z})$.

Based on the above two relaxations, the problem $\mathrm{P}3$ with $\widetilde{\mathbf{v}}\in\mathcal{F}_{(z_{1},\cdots,z_{N})}$ is relaxed into

[TABLE]

Using the result from (17)-(18), the optimal $\{s_{k}\}$ to the problem (52) is given by

[TABLE]

with $\delta$ obtained from (24). Putting (53) into the objective function of (52), we immediately obtain $\Psi$ in (23), which is obviously a lower bound to the objective function $\Xi\left(\widetilde{\mathbf{v}}\right)$ of $\mathrm{P}3$ for any $\widetilde{\mathbf{v}}\in\mathcal{F}_{(z_{1},\cdots,z_{N})}$.