# Game-Theoretic Algorithms for Conditional Moment Matching

Gokul Swamy  Sanjiban Choudhury  J. Andrew Bagnell  Zhiwei Steven Wu

###### Abstract

A variety of problems in econometrics and machine learning, including
instrumental variable regression and Bellman residual minimization, can be
formulated as satisfying a set of conditional moment restrictions (CMR). We
derive a general, game-theoretic strategy for satisfying CMR that scales to
nonlinear problems, is amenable to gradient-based optimization, and is able to
account for finite sample uncertainty. We recover the approaches of (Dikkala
et al., 2020) and (Dai et al., 2018) as special cases of our general framework
before detailing various extensions and how to efficiently solve the game
defined by CMR.

Machine Learning, ICML

##  1 Introduction

Let $X$, $Y$, and $Z$ be random variables on (potentially non-finite) sample
spaces $\mathcal{X}$, $\mathcal{Y}$, and $\mathcal{Z}$. We are interested in
finding a function $h\in\mathcal{H}$ that satisfies a set of conditional
moment restrictions (Chamberlain, 1987),

$\forall z\in\mathcal{Z},\,\mathbb{E}[Y|z]=\mathbb{E}[h(X)|z],$ (1)

or CMR for short. While this problem might seem a bit abstract, it is at the
core of two disparate problems in machine learning: instrumental variable
regression and Bellman-residual minimization. We give a brief introduction to
each before presenting a unified method for solving for consistent $h$ from
finite samples that elegantly scales to nonlinear problems. Throughout, we
assume we optimize over a class $\mathcal{H}$ that is convex, compact, closed
under negation, and of finite Rademacher complexity.

Our key insight is that we can formulate conditional moment matching as a
zero-sum game, allowing us to both eliminate double sample issues and
explicitly reason about the effects of constraint relaxation. We call this
family of techniques estimation via relaxation.

###  1.1 Instrumental Variable Regression

Let’s assume that $X$, $Y$, and $Z$ have the following dependency structure:

$Z$$X$$Y$$U$$g$$h$ Figure 1: The graphical model considered in instrumental
variable regression. We are interested in finding $h$, the causal relationship
from $X$ to $Y$, even though there is an unobserved confounder, $U$.

We refer to $X$ as the treatment and $Y$ as the response or outcome. Given a
dataset of $(x,y,z)$ tuples, we are interested in determining the causal
relationship between $X$ and $Y$, $\mathbb{E}[Y|\text{do}(x)]$. In the above
graphical model, this is equivalent to determining $h$. Because of the
presence of an unobserved confounder, $U$, that affects both $X$ and $Y$,
standard regression (e.g. Ordinary Least Squares or OLS) can produce
inconsistent estimates. If we only have observational data and are unable to
perform randomized control trials, a canonical technique to recover $h$ is
Instrumental Variable Regression (IVR) (Winship & Morgan, 1999). Formally, an
instrument $Z$ must satisfy:

  1. 1.

Unconfounded Instrument: $Z\perp\\!\\!\\!\perp U$ – i.e. “independent”
randomization from instrument.

  2. 2.

Exclusion: $Z\perp\\!\\!\\!\perp Y|X$ – i.e. no extraneous paths.

  3. 3.

Relevance: $Z\not\\!\perp\\!\\!\\!\perp X$ – i.e. conditioning has an effect.

We note that the unconfounded instrument and exclusion conditions are
structural assumptions that cannot be checked from observational data.
Assuming access to such a $Z$ and linear relationships between all variables,
one can recover $h(x)=\beta x$ by computing
$\beta=\mathbb{E}[ZX]^{-1}\mathbb{E}[ZY]$. Notice that the recovered $h$
satisfies the CMR (1). With a bit of linear algebra, one can show this
calculation is equivalent to a Two-Stage Ordinary Least Squares (2SLS)
procedure, for which one first regresses from $Z$ to $X$ and then regresses
from the predicted $\hat{X}$ to $Y$, returning the latter coefficients.

We focus on the general, nonlinear problem, which can still be formulated in
terms of CMRs. We assume that noise $U$ enters additively to $Y$ 111This is an
assumption inherited from standard regression. Without it, one can only bound
the treatment effect (Kilbertus et al., 2020)., and write out the following
equations:

$X=g(Z,U),\vspace{-4pt}$ (2) $Y=h(X)+U.\vspace{-4pt}$ (3)

Without loss of generality, we assume that $\mathbb{E}[U]=0$. As is standard
in the causal inference literature (Nie et al., 2020), we assume that our
collected data satisfies consistency of potential outcomes (i.e. the data
generating process did not change since our dataset was collected (Rehkopf et
al., 2016)) and overlap (i.e. we have data for all $X$ we care about (Imbens,
2000)). Perhaps the most natural way to tackle this problem would be to
minimize some notion of squared error between the two sides of the CMR:

$\min_{h\in\mathcal{H}}\mathbb{E}_{z}[\mathbb{E}[Y-h(x)|z]^{2}].\vspace{-4pt}$
(4)

For nonlinear problems, gradient-based optimization is a common and scalable
technique. Unfortunately, differentiating the squared expectation in the
preceding expression leads to a “double sample” issue in which one needs two
independent samples from $P(X|z)$ to compute each gradient step, which can be
quite challenging outside of simulated environments. To resolve this issue,
(Dikkala et al., 2020) propose instead solving the following zero-sum game:

$\min_{h\in\mathcal{H}}\max_{f}\mathbb{E}[2(Y-h(X))f(Z)-f(Z)^{2}],$ (5)

which they prove produces consistent estimates of $h$ and does not involve a
squared expectation. We show in Sec. 1 that this is an example of a general
technique for satisfying CMR.

###  1.2 Bellman Residual Minimization

Consider a Markov Decision Process (MDP) parameterized by
$\langle\mathcal{S},\mathcal{A},\mathcal{T},r,\gamma\rangle$, where
$\mathcal{S}$ is the state space, $\mathcal{A}$ is the action space,
$\mathcal{T}:\mathcal{S}\times\mathcal{A}\rightarrow\Delta(\mathcal{S})$ is
the transition operator, $r:\mathcal{S}\times\mathcal{A}\rightarrow[-1,1]$ is
a reward function, and $\gamma$ is the discount factor. It is a well-known
fact that one can find the optimal policy, $\pi$, for this MDP by first
computing the optimal value function via the Bellman Equation,

$V(s)=\max_{a\in\mathcal{A}}r(s,a)+\gamma\mathbb{E}_{s^{\prime}\sim\mathcal{T}(s,a)}[V(s^{\prime})],$
(6)

and choosing actions greedily (Puterman, 2014):

$\pi(s)=\arg\max_{a\in\mathcal{A}}r(s,a)+\gamma\mathbb{E}_{s^{\prime}\sim\mathcal{T}(s,a)}[V(s^{\prime})].$
(7)

We can combine these two expressions into one that is satisfied by the optimal
$(V,\pi)$ pair:

$V(s)=\mathbb{E}_{\begin{subarray}{c}a\sim\pi(s)\\\
s^{\prime}\sim\mathcal{T}(s,a)\end{subarray}}[r(s,a)+\gamma V(s^{\prime})].$
(8)

When one cannot simply enumerate all states and actions to perform either
policy or value iteration (Russell & Norvig, 2002), a standard technique is
Bellman residual minimization – minimizing the difference between the two
sides of (8). The updated value function can then be used to compute a more
optimal policy via the greedy optimization of (7). We can write (8) in the
template of CMR by setting

$Z=s,a\vspace{-4pt}$ (9) $X=s,\,h(X)=V(s)\vspace{-4pt}$ (10)
$Y=r(s,a)+\gamma\mathbb{E}_{s^{\prime}\sim\mathcal{T}(s,a)}[V(s^{\prime})].\vspace{-4pt}$
(11)

Similar to (4), one could attempt to satisfy the CMR by minimizing the
expected squared Bellman error:

$\min_{V:\mathcal{S}\rightarrow\mathbb{R}}\mathbb{E}[\mathbb{E}[V(s)-(r(s,a)+\gamma\mathbb{E}_{s^{\prime}\sim\mathcal{T}(s,a)}[V(s^{\prime})])|s,a]^{2}],\vspace{-4pt}$
(12)

which also has double sample issues (Baird, 1995). Recently, (Dai et al.,
2018) suggested that one could avoid double sample issues by instead solving
the following zero-sum game:222For simplicity, we drop the entropy
regularization in their objective. Our derivations can be easily extended.

$\min_{V:\mathcal{S}\rightarrow\mathbb{R}}\max_{f:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}}\mathbb{E}_{s,a,s^{\prime}}[2(r(s,a)+V(s^{\prime})-V(s))f(s,a)-f(s,a)^{2}],\vspace{-4pt}$
(13)

which they derive via an appeal to convex conjugates (Boyd et al., 2004). We
present an alternative construction of this objective which helps elucidate
the properties of the algorithm proposed by (Dai et al., 2018).

##  2 Estimation via Relaxation

Setting | $Z$ |  $X$ and $h(X)$ | $Y$
---|---|---|---
IVR | Instrument | Treatment | Outcome
RL | $(s,a)$ |  $s$ and $V(s)$ | $r(s,a)+\gamma\mathbb{E}_{s^{\prime}}[V(s^{\prime})]$
Table 1: The work of both (Dikkala et al., 2020) and (Dai et al., 2018) can be
seen as examples of a more general template for satisfying conditional moment
restrictions from finite samples and for nonlinear problems.

Throughout this section, we use bold fonts to designate vectors. To aid us in
our quest to satisfy the conditional moment restrictions, we are given access
to $N$ samples from the joint distribution of $X$, $Y$, and $Z$. Because we
have finite samples and can therefore only estimate conditional expectations
up to some tolerance, it is natural to relax the CMRs to

$\begin{array}[]{ll@{}ll}\min_{h\in\mathcal{H},\,\boldsymbol{\delta}}&J(\boldsymbol{\delta})\\\
\text{s.t.}&|\mathop{{}\mathbb{E}}[Y-h(X)|z]|\leq\delta_{z}\quad\forall
z\in\mathcal{Z},\\\ \end{array}$ (14)

where the $\delta_{z}$ are slack variables and $J$ is some convex function of
$\boldsymbol{\delta}$ that keeps the slacks from getting too large. We note
that (14) is reminiscent of the classical framework of Tikhinov regularization
before presenting some interesting properties that result from setting
$J(\boldsymbol{\delta})=\frac{1}{2}\mathbb{E}_{z}[\delta_{z}^{2}]$. The
Lagrangian with a $P(z)$-weighted inner product is

$L(h,\boldsymbol{\delta},\boldsymbol{\lambda})=\sum_{z\in\mathcal{Z}}P(z)\lambda_{z}(\mathbb{E}[Y-h(X)|z]-\delta_{z})+P(z)\frac{1}{2}\delta_{z}^{2},$
(15)

where $\boldsymbol{\lambda}$ are the Lagrange multipliers.333While written
above in terms of finite $z$, our derivation easily extend to infinite sets.
Applying the stationarity component of the KKT conditions, one arrives at

$\nabla_{\delta_{z}}L(h,\boldsymbol{\delta},\boldsymbol{\lambda})=-P(z)\lambda_{z}+P(z)\delta_{z}=0,$
(16)

implying that $\delta_{z}=\lambda_{z}$. Plugging this back into the
Lagrangian, we can simplify our function to

$L(h,\boldsymbol{\lambda})=\sum_{z\in\mathcal{Z}}P(z)\lambda_{z}\mathbb{E}[Y-h(X)|z]-P(z)\frac{1}{2}\lambda_{z}^{2}.$
(17)

We refer to (17) as the Regularized Lagrangian or ReLa for short. Now, solving
for the optimal Lagrange multipliers via stationarity, we arrive at

$\nabla_{\lambda_{z}}L(h,\boldsymbol{\lambda})=P(z)\mathbb{E}[Y-h(X)|z]-P(z)\lambda_{z}=0,$
(18)

which implies the equilibrium $\lambda_{z}$ is equal to
$\mathbb{E}[Y-h(X)|z]$. Plugging this back into (17) recovers function

$L(h)=\sum_{z\in\mathcal{Z}}P(z)\mathbb{E}[Y-h(X)|z]^{2}.$ (19)

Thus, in the population limit, we are optimizing the conditional MSE of (Chen
& Pouzo, 2012), leading to consistent estimates.

###  2.1 Generative Modeling Approach

Perhaps the most immediate way to minimize (19) over $h\in\mathcal{H}$ would
be to minimize the empirical MSE,

$\min_{h\in\mathcal{H}}\frac{1}{N}\sum_{i}^{N}(y_{i}-h(x_{i}))^{2}.$ (20)

Unfortunately, this only gives us a function that matches unconditional
moments (i.e. $\mathbb{E}[Y]=\mathbb{E}[h(X)]$). For the IVR setting, this
would give us the inconsistent, OLS-like estimates of $h$ we are explicitly
trying to avoid. For the MDP setting, this would not produce a valid value
function.

Instead, one could learn the distribution $P(X|z)=g(z)$ and pass samples from
it to a candidate $h$, ensuring one is attempting to match the conditional
moments. Intuitively, this is the generalization of the 2SLS procedure to
nonlinear functions. Because the second stage is nonlinear, one cannot simply
compute the first moment of the $P(X|z)$ distribution (which is recovered by
linearly regressing from $X$ to $Z$ in the 2SLS procedure). To see this,
consider $\mathcal{H}$ being the set of quadratic functions of $X$. Then,

$\displaystyle\mathbb{E}[h(X)|z]$
$\displaystyle=\sum_{x\in\mathcal{X}}p(x|z)h(x)$ (21)
$\displaystyle=\sum_{x\in\mathcal{X}}p(x|z)(ax^{2}+bx+c)$ (22)
$\displaystyle=a\mathbb{E}[X^{2}|z]+b\mathbb{E}[X|z]+c.$ (23)

Here, one needs to have access to $\mathbb{E}[X^{2}|z]$ to check the CMR. For
more complex $\mathcal{H}$, one therefore needs to learn the entire $P(X|z)$.
This kind of approach was first proposed for the IVR setting by (Hartford et
al., 2017) and amounts to first learning a $g(z)$ via maximum likelihood
estimation and then solving

$\min_{h\in\mathcal{H}}\frac{1}{N}\sum_{i}^{N}(y_{i}-\mathbb{E}_{\hat{x}\sim
g(z_{i})}[h(\hat{x})])^{2}.$ (24)

Unfortunately, this approach suffers from the well-known “double-sample” issue
where multiple independent samples from $g(z)$ are required to compute
gradients of $h$. To see this, note that the gradient w.r.t. $h$ of (24) is

$\sum_{i}^{N}(y_{i}-\mathbb{E}_{\hat{x}\sim\hat{g}(z_{i})}[h(\hat{x})])(-\mathbb{E}_{\hat{x}\sim\hat{g}(z_{i})}[\frac{\partial}{\partial
h}h(\hat{x})]).$ (25)

Approximating this gradient with a single sample,

$\sum_{i}^{N}(\mathop{{}\mathbb{E}}_{\hat{x}\sim\hat{g}(z_{i})}[(y_{i}-h(\hat{x}))\frac{\partial}{\partial
h}h(\hat{x})]),$ (26)

can produce biased estimates of the gradient that are inconsistent in the
limit of infinite data. Additionally, learning an accurate first-stage model
might be quite challenging for some problems.

###  2.2 Game-Theoretic Approach

Ideally, we would like to avoid the added complexity of learning a generative
model and the double-sampling required for gradient based-optimization. One
gets a two-for-one deal by instead solving the two-player zero-sum game with
the ReLa (17) as the payoff. Denoting by
$f\in\mathcal{F}\subseteq\\{\mathcal{Z}\rightarrow\mathbb{R}\\}$ the function
that maps $z$’s to corresponding Lagrange multipliers, we can write this game
as:

$\min_{h\in\mathcal{H}}\max_{f\in\mathcal{F}}\mathbb{E}[2(Y-h(X))f(Z)-f(Z)^{2}].$
(27)

Notice the similarity of this expression to the objectives of (Dikkala et al.,
2020) (5) and (Dai et al., 2018) (13) and that there are no squared
expectations of $h$ and therefore no double sample issues. Additionally, one
does not need to learn a generative model of $P(X|z)$ for these sorts of game-
theoretic approaches. Under the assumption that
$\mathbb{E}[Y-h(X)|z]\in\mathcal{F}$, finding a Nash equilibrium of this game
corresponds to finding a $h\in\mathcal{H}$ that is close in an $\ell_{2}$
sense to satisfying the CMR.

###  2.3 Takeaways Thus Far

We pause briefly to consider some of the important facts revealed to us by the
above derivation:

  1. 1.

The approaches of (Dikkala et al., 2020) and (Dai et al., 2018) are solving a
relaxed CMR problem with a penalty on the $\ell_{2}$ norm of the constraint
violation. Thus, their consistency is in an $\ell_{2}$ sense rather than a
uniform, $\ell_{\infty}$ sense.

  2. 2.

The choice of $J$ allows us to control how slack is dispersed between the CMRs
– derivations like the above allow us to explicitly consider the effects of
the kind of regularization we choose.

  3. 3.

The value of $f(z)$ (regardless of the choice of $J$) is exactly the amount of
conditional moment mismatch the optimization procedure tolerates. If this
value is particularly high for values of $z$ where we have many samples, we
might need to consider a more expressive $\mathcal{H}$. Thus, we can
explicitly reason about the effects of finite-sample uncertainty in satisfying
CMRs.

##  3 Extensions

We now briefly sketch two extensions to our above setup before diving into how
to solve ReLa games efficiently.

###  3.1 Constraints on $\boldsymbol{\delta}$

A natural question after the above discussion might be how one could include
information about sample uncertainty into the optimization procedure, rather
than performing an after-the-fact check. Let $n_{z}$ denote the number of
samples we have for $z$. For simplicity, we assume $n_{z}>0,\,\forall
z\in\mathcal{Z}$. A Hoeffding bound tells us that we should expect a sample
expectation $\hat{\mathbb{E}}_{n_{z}}[\cdot|z]$ to be within
$\propto\frac{1}{\sqrt{n_{z}}}$ of the population expectation. We refer to
this sampling error as $\epsilon_{z}$. Then, in expectation over $Z$,

$\mathbb{E}_{Z}[\epsilon_{z}^{2}]=\sum_{z\in\mathcal{Z}}P(z)\epsilon_{z}^{2}\approx\sum_{z\in\mathcal{Z}}\frac{n_{z}}{N}(\frac{1}{\sqrt{n_{z}}})^{2}=\frac{|\mathcal{Z}|}{N}=\kappa(N).$
(28)

Our relaxed optimization problem then becomes

$\begin{array}[]{ll@{}ll}\min_{h\in\mathcal{H},\,\boldsymbol{\delta}}&0\\\
\text{s.t.}&|\mathop{{}\mathbb{E}}[Y-h(X)|z]|\leq\delta_{z}\quad\forall
z\in\mathcal{Z}\\\ &\mathbb{E}_{Z}[\delta_{z}^{2}]\leq\kappa(N).\end{array}$
(29)

This form, reminiscent of Ivanov regularization, takes into account sample
uncertainty explicitly and could be solved via standard constrained
optimization machinery like Augmented Lagrangians (Hestenes, 1969). We note
that in practice, it is far more common to solve the Tikhanov-style problem
(14) and scale $J$ based on performance on some holdout data.

###  3.2 Regularization of $h$

Consider two solutions $(h_{1},f_{1})$, $(h_{2},f_{2})$ of the ReLa Game (27)
such that

$\displaystyle|\mathbb{E}[Y-h_{1}(X)|z]|\leq f_{1}(z),\,\forall
z\in\mathcal{Z},$ (30) $\displaystyle|\mathbb{E}[Y-h_{2}(X)|z]|\leq
f_{2}(z),\,\forall z\in\mathcal{Z},$ (31)
$\displaystyle\mathbb{E}_{Z}[f_{1}(z)^{2}]=\mathbb{E}_{Z}[f_{2}(z)^{2}].$ (32)

Our game-theoretic perspective wouldn’t be able to break ties between $h_{1}$
and $h_{2}$, regardless of desirable properties $h_{1}$ might have over
$h_{2}$ like smoothness in $X$ or being the max-margin classifier. This is
because as long as two solutions exist within CMR-violation balls of the same
size (e.g. in $\ell_{2}$ norm), we consider them equally optimal. Thus, some
sort of regularization on $h$, producing an optimization problem of the form

$\begin{array}[]{ll@{}ll}\min_{h\in\mathcal{H},\,\boldsymbol{\delta}}&J(\boldsymbol{\delta})+\alpha
R(h)\\\
\text{s.t.}&|\mathop{{}\mathbb{E}}[Y-h(X)|z]|\leq\delta_{z}\quad\forall
z\in\mathcal{Z}\\\ \end{array}$ (33)

with $\alpha$ tuned empirically on hold-out data might be helpful for finding
$h$ with additional desirable properties. Concretely, one could regularize to
the OLS solution by setting $R(h)=\mathbb{E}[(h(X)-Y)^{2}]$, encouraging $h$
to also match unconditional moments. Both (Dikkala et al., 2020) and (Dai et
al., 2018) consider regularization on $h$ (e.g. in
$||\cdot||^{2}_{\mathcal{H}}$ norm or by entropy-regularizing the policy),
which we identify as a form of tie-breaking.

###  3.3 Efficiently Solving the ReLa Game

Input: Dataset $\mathcal{D}$ of $(x,y,z)$ tuples. No-regret algorithm over
$\mathcal{H}$, Best-response oracle over $\mathcal{F}$, Threshold $\epsilon$

Output: Causal effect of $X$ on $Y$, $h$

Set $t=1$, $h^{t}\in\mathcal{H}$, $f^{t}\in\mathcal{F}$,
$L(h^{t},f^{t})=2\epsilon$

while $L(h^{t},f^{t})>\epsilon$ do

$L(h^{t},f^{t})=\mathop{{}\mathbb{E}}_{(x,y,z)\sim\mathcal{D}}[2(y-h^{t}(x))f^{t}(z)-f^{t}(z)^{2}]$

No-regret alg. computes $h^{t}$ over $L(\cdot,f^{t})$ history.

Best-response computes $f^{t}=\arg\max_{f\in\mathcal{F}}L(h^{t},\cdot)$.

$t\leftarrow t+1$

end while

Return $h^{t}$.

Algorithm 1 No-Regret Conditional Moment Matching

Solving a two-player game like (27) can be done provably efficiently via a
reduction to no-regret online learning, following the classic analysis of
(Freund & Schapire, 1997). For completeness, we provide such a procedure in
Algorithm 1 and prove efficiency below.

######  Proof.

Let $L(h,f)=\mathbb{E}[2(Y-h(X))f(Z)-f(Z)^{2}]$. Our no-regret assumption
tells us that

$\frac{1}{N}\sum_{t}^{N}L(h^{t},f^{t})-\frac{1}{N}\min_{h\in\mathcal{H}}\sum_{t}^{N}L(h,f^{t})\leq\epsilon.$
(34)

for some $N$ that is $\text{poly}(\frac{1}{\epsilon})$. Under a realizability
assumption, we can write that

$\frac{1}{N}\sum_{t}^{N}L(h^{t},f^{t})\leq\epsilon+\frac{1}{N}\min_{h\in\mathcal{H}}\sum_{t}^{N}L(h,f^{t})\leq\epsilon.$
(35)

Then, utilizing the fact that there must be at least one element in an average
that is at most the value of the average,

$\min_{t}L(h^{t},f^{t})\leq\frac{1}{N}\sum_{t}^{N}L(h^{t},f^{t})\leq\epsilon.$
(36)

To complete the proof, we recall that $f^{t}$ is chosen as the best response
to $h^{t}$ in Algorithm 1, giving us that:

$\min_{t}\max_{f}L(h^{t},f)\leq\epsilon$ (37)

Thus, the $h^{t}$ that minimizes $L(h^{t},f^{t})$ is half of an
$\epsilon$-approximate Nash equilibrium of the ReLa game and can be computed
within $\text{poly}(\frac{1}{\epsilon})$ iterations. In our setting,
$\epsilon$ corresponds to the additional expected CMR violation our recovered
$h$ suffers on top of the best element in $\mathcal{H}$.

We note this reduction allows one to plug in any no-regret algorithm (e.g.
Follow the Regularized Leader (McMahan, 2011) or Multiplicative Weights (Arora
et al., 2012)) and the efficiency and approximate CMR satisfaction result to
still hold. Practically, one could instantiate this no-regret reduction via a
GAN-like optimization procedure with the learning rate for the $h$ player much
lower than that of the $f$ player, simulating the no-regret vs. best response
iterations of Algorithm 1. While the connection to no-regret online learning
is not made explicit, both (Dikkala et al., 2020) and (Dai et al., 2018)
follow an approach in this vein to scale their methods to high-dimensional
tasks.

##  4 Discussion

We unify the techniques of (Dikkala et al., 2020) and (Dai et al., 2018) under
the umbrella of satisfying relaxed conditional moment restrictions, providing
a framework that allows one to reason about sample uncertainty and the effects
of regularization on such problems. We further consider explicit constraints
on slack variables, the tie-breaking benefits of regularizing $h$, and provide
a reduction to no-regret online learning for efficient solving of Regularized
Lagrangian Games.

Moving forward, we would be interested in passing the higher bar of matching
not just the first moments of $Y|z$ but the whole distribution, as well as
fitting other problems into the framework of conditional moment restrictions.

##  5 Acknowledgements

We thank Keegan Harris and Logan Stapleton for comments on a draft of this
work. ZSW was supported in part by the NSF FAI Award #1939606, a Google
Faculty Research Award, a J.P. Morgan Faculty Award, a Facebook Research
Award, and a Mozilla Research Grant. GS is supported by NSF Award #1952085 and
his family and friends.

## References

  * Arora et al. (2012) Arora, S., Hazan, E., and Kale, S.  The multiplicative weights update method: a meta-algorithm and applications.  _Theory of Computing_ , 8(1):121–164, 2012.
  * Baird (1995) Baird, L.  Residual algorithms: Reinforcement learning with function approximation.  In _Machine Learning Proceedings 1995_ , pp. 30–37. Elsevier, 1995\.
  * Boyd et al. (2004) Boyd, S., Boyd, S. P., and Vandenberghe, L.  _Convex optimization_.  Cambridge university press, 2004.
  * Chamberlain (1987) Chamberlain, G.  Asymptotic efficiency in estimation with conditional moment restrictions.  _Journal of econometrics_ , 34(3):305–334, 1987\.
  * Chen & Pouzo (2012) Chen, X. and Pouzo, D.  Estimation of nonparametric conditional moment models with possibly nonsmooth generalized residuals.  _Econometrica_ , 80(1):277–321, 2012.
  * Dai et al. (2018) Dai, B., Shaw, A., Li, L., Xiao, L., He, N., Liu, Z., Chen, J., and Song, L.  Sbeed: Convergent reinforcement learning with nonlinear function approximation, 2018.
  * Dikkala et al. (2020) Dikkala, N., Lewis, G., Mackey, L., and Syrgkanis, V.  Minimax estimation of conditional moment models, 2020.
  * Freund & Schapire (1997) Freund, Y. and Schapire, R. E.  A decision-theoretic generalization of on-line learning and an application to boosting.  _Journal of computer and system sciences_ , 55(1):119–139, 1997.
  * Hartford et al. (2017) Hartford, J., Lewis, G., Leyton-Brown, K., and Taddy, M.  Deep iv: A flexible approach for counterfactual prediction.  In _International Conference on Machine Learning_ , pp. 1414–1423. PMLR, 2017.
  * Hestenes (1969) Hestenes, M. R.  Multiplier and gradient methods.  _Journal of optimization theory and applications_ , 4(5):303–320, 1969.
  * Imbens (2000) Imbens, G.  The role of the propensity score in estimating dose-response functions.  _Biometrika_ , 87(3):706–710, 2000.
  * Kilbertus et al. (2020) Kilbertus, N., Kusner, M. J., and Silva, R.  A class of algorithms for general instrumental variable models, 2020.
  * McMahan (2011) McMahan, H. B.  Follow-the-regularized-leader and mirror descent: Equivalence theorems and l1 regularization.  In _Proceedings of the 14th International Conference on Artificial Intelligence and Statistics (AISTATS)_ , 2011.
  * Nie et al. (2020) Nie, X., Brunskill, E., and Wager, S.  Learning when-to-treat policies, 2020.
  * Puterman (2014) Puterman, M. L.  _Markov decision processes: discrete stochastic dynamic programming_.  John Wiley & Sons, 2014.
  * Rehkopf et al. (2016) Rehkopf, D. H., Glymour, M. M., and Osypuk, T. L.  The consistency assumption for causal inference in social epidemiology: when a rose is not a rose.  _Current epidemiology reports_ , 3(1):63–71, 2016\.
  * Russell & Norvig (2002) Russell, S. and Norvig, P.  Artificial intelligence: a modern approach.  2002\.
  * Winship & Morgan (1999) Winship, C. and Morgan, S. L.  The estimation of causal effects from observational data.  _Annual review of sociology_ , 25(1):659–706, 1999.


# On Energy Allocation and Data Scheduling in Backscatter Networks with Multi-
antenna Readers

Mohammad Movahednasab, Mohammad Reza Pakravan, Behrooz Makki, and Tommy
Svensson M. Movahednasab and M. R. Pakravan are with Sharif University of
Technology, Tehran, Iran (Emails<EMAIL_ADDRESS>pakravan@sharif.edu). Behrooz Makki is in Gothenburg, Sweden (Email:
behrooz.makki@gmail.com). Tommy Svensson is with Chalmers University of
Technology, Gothenburg, Sweden (Email: tommy.svensson@chalmers.se).

###### Abstract

In this paper, we study the throughput utility functions in buffer-equipped
monostatic backscatter communication networks with multi-antenna Readers. In
the considered model, the backscatter nodes (BNs) store the data in their
buffers before transmission to the Reader. We investigate three utility
functions, namely, the sum, the proportional and the common throughput. We
design online admission policies, corresponding to each utility function, to
determine how much data can be admitted in the buffers. Moreover, we propose
an online data link control policy for jointly controlling the transmit and
receive beamforming vectors as well as the reflection coefficients of the BNs.
The proposed policies for data admission and data link control jointly
optimize the throughput utility, while stabilizing the buffers. We adopt the
min-drift-plus-penalty (MDPP) method in designing the control policies.
Following the MDPP method, we cast the optimal data link control and the data
admission policies as solutions of two independent optimization problems which
should be solved in each time slot. The optimization problem corresponding to
the data link control is non-convex and does not have a trivial solution.
Using Lagrangian dual and quadratic transforms, we find a closed-form
iterative solution. Finally, we use the results on the achievable rates of
finite blocklength codes to study the system performance in the cases with
short packets. As demonstrated, the proposed policies achieve optimal utility
and stabilize the data buffers in the BNs.

###### Index Terms:

Backscatter communication, radio frequency identification, fairness, min-
drift-plus-penalty, Lyapunov optimization, max-min throughput, proportional
throughput, sum throughput, finite blocklength analysis, wireless energy
transfer, energy harvesting, green communications, Internet of things, IoT.

##  I INTRODUCTION

With the emergence of the Internet of things (IoT) era, the number of
connected devices is increasing rapidly. It is predicted that there will be
around 5 billion connected devices in 2025 [1], a great number of which will
be portable and low-power. This explosion of the low-powered devices calls for
replacing batteries with new energy sources to ensure continuous operation of
devices. The major challenges with the battery-powered devices are the
increase of the devices’ form factor and the high cost for recharging and
replacement of the batteries [2]. Moreover, in some applications such as
biomedical implants inside human bodies or distributed monitoring sensors,
replacing the batteries may be infeasible [3, 4, 5].

Backscatter communication networks (BCNs) are considered to be a prominent
solution for low-power and low-cost communications. A BCN compromises a Reader
and, possibly, multiple backscattering nodes (BNs) with most bulky and active
communication modules moved to the Reader. The BNs transmit data to the Reader
via reflecting and modulating the incident radio frequency signal by adapting
the level of antenna mismatch to vary the reflection coefficient [4]. Based on
the source of the radio frequency signal, which supplies the required energy
for communication, three configurations for the BCN, namely, monostatic,
bistatic and ambient, are considered. In monostatic configuration, the Reader
emits a carrier and receives the backscattered data, while in bistatic
configuration one or several power beacons emit carries rather than the Reader
itself [6]. Moreover, in ambient BCNs there is no dedicated energy transmitter
and the BNs backscatter the existing radio frequency signals, e.g., WiFi or
digital television signals [7].

The low-energy transmission efficiency is a major problem in BCNs. However,
the studies in, e.g., [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19], show
that exploiting multiple antennas increases the energy efficiency remarkably.
Considering the energy required for channel training, [8] optimizes the
transmit beamforming to maximize the harvested energy by the BNs. In [9], a
blind adaptive beamforming scheme is introduced to increase the reading range
of the radio frequency identification tags. Also, [10] and [11] propose low
complexity algorithms for optimizing the transmit and receive beamforming to
maximize the sum and minimum throughput of all BNs, respectively. Considering
ambient BCNs, [12, 13, 14, 15] study the receive beamforming optimization of a
multi-antenna Reader. Whereas, [16] and [17] study optimal detectors for
ambient BCNs with multi-antenna BNs and a single antenna Reader. In [18] and
[19], the achievable diversity order of a multiple-input-multiple-output
monostatic BCN is studied. Furthermore, optimizing the reflection coefficients
of the BNs is studied in, e.g., [20, 21, 22, 23, 24, 25], where maximizing
energy efficiency, throughput or fairness in the BCNs is investigated.

Along with energy efficiency, one of the challenges of the BCNs is the
unpredictability of the channel state and the available energy in ambient
configuration, which makes the optimal scheduling difficult. To tackle this
problem, stochastic approaches are adopted in [26, 27, 28, 29, 30] to design
online control algorithms. Specifically, in [26, 27, 28], long-term throughput
optimization in ambient BCNs is studied through reinforcement learning
methods. Whereas, in [29] the authors use reinforcement learning to propose a
data admission and data scheduling policy for a monostatic BCN. Finally, [30]
uses min-drift-plus-penalty (MDPP) method to maximize throughput in an ambient
BCN.

In this work, we concentrate on optimizing different throughput utility
functions in monostatic BCNs with multi-antenna Readers. In our considered
model, the BNs store the data in their buffers before transmission to the
Reader. We investigate three different utility functions including the sum,
the proportional and the common throughput. We design three online admission
policies, corresponding to each throughput utility function, to determine how
much data can be admitted in the buffers in each time slot. Moreover, we
propose an online data link control policy for jointly controlling the
transmit and receive beamforming vectors as well as the reflection
coefficients of the BNs. The proposed policies for data admission and data
link control jointly optimize the throughput utility, while stabilizing the
buffers.

Considering the channel state randomness, we adopt the MDPP method in
designing the control policies. Following the MDPP method, we cast the optimal
data link control and the data admission policies as solutions of two
independent optimization problems which should be solved in each time slot.
The optimization problem corresponding to the data link control is non-convex
and does not have a trivial solution. We transform this problem into an
equivalent form, which makes it possible to find a closed-form iterative
algorithm that improves the utility in each iteration. Furthermore,
considering each utility function, we solve the corresponding data admission
problem and find closed-form solutions. Finally, we use the results on the
achievable rates of finite blocklength codes [31, 32, 33, 34] to study the
system performance in the cases with short packets.

The differences in the considered system model and problem formulation makes
our paper different from those in the literature. Specifically, this paper is
different from [12, 13, 14, 15, 16, 17, 20, 21, 22, 24, 26, 27, 28, 29, 30],
because we study a monostatic BCN, where we jointly optimize the transmit and
receive beamforming vectors. As opposed to [8, 9, 10, 11, 12, 13, 14, 15, 16,
17, 18, 19, 20, 21, 22, 23, 24, 25], we consider BNs equipped with buffers,
and hence follow an stochastic approach to propose an online policy to
maximize the throughput utilities and stabilize the buffers. Moreover,
different from the state-of-the-art works, we study different throughput
utilities in a unified framewok, perform finite block-length analysis, and
compare the results in terms of fairness.

Our analytical and simulations results show the significance of our proposed
control policy. Specifically, we show in Theorem LABEL:th:optimal that, under
the optimal solutions of the link scheduling and data admission problems, the
average throughput utility will be within $\mathcal{O}({1\over V})$ of the
optimal utility, where $V$ is a control parameter. Moreover, the level of the
stored data in the buffers will be kept under a certain level of
$\mathcal{O}({V})$. Furthermore, we compare our proposed scheme with the
benchmarks [10] and [11], which verifies the optimality of our proposed
policy. Particularly, [10] and [11] propose optimal control policies to
maximize the sum and the minimum channel rate among BNs in a time slot base,
respectively. However, the BNs in these two papers do not consider buffers. We
show that the sum throughput utility under our proposed policy achieves the
maximum sum channel rate in [10], while stabilizing the data buffers. Whereas,
our proposed policy for optimizing the common throughput improves the result
in [11], which shows the significant effect of adopting buffers in the BNs.

The rest of the paper is organized as follows. The considered system model and
our problem formulation are illustrated in Section II. The proposed control
policy as well as its performance analysis are presented in Section III.
Simulation results are presented in Section IV, and finally, Section V
concludes the paper.

Notation: Matrices and vectors are denoted by small and capital boldface
letters, respectively. Moreover, unless otherwise mentioned, vectors are
single-column matrices. Also, $(.)^{T}$ , $(.)^{H}$ and $(.)^{\ast}$ denote
transpose, conjugate transpose and elementwise conjugate of a matrix,
respectively. Then, $|.|$ denotes the absolute value (or the modulus for
complex numbers) and $\lVert.\lVert$ denotes the norm of vectors. Finally,
$\mathbb{R}$ and $\mathbb{C}$ represent the real and complex number sets,
respectively, $\operatorname{Re}\\{.\\}$ is used to denote the real part of a
complex number, $\operatorname{E}\\{.\\}$ is the expectation, $I_{N}$
represents the $N\times N$ identity matrix, and
$[a]^{+}=\max\\{a,0\\},\;\forall a\in\mathbb{R}$.

##  II System Model and Problem Formulation

An example of the considered network structure is depicted in Fig. 1. We
consider a BCN consisting of a multi-antenna Reader and $N$ single-antenna
BNs. Let $\text{BN}_{n},\;\forall n\in\mathcal{N}\triangleq\\{1,\ldots,N\\}$,
denote the $n$-th BN in the network. The BNs transmit the data stored in their
buffers to the Reader through backscattering the energy emitted by the Reader.
While the BNs rely on the Reader energy for data transmission, they have an
internal battery which powers up their low-power circuits. Accordingly, the
BNs are semi-passive devices and are able to handle data sensing or other
internal tasks without the Reader energy [35, 36, 37]. The Reader is equipped
with $M$ antennas to focus its emitted energy towards the intended BNs.
Moreover, the presence of multiple antennas enables the Reader to receive data
from multiple BNs through receive beamforming technique. The Reader adopts a
kind of full duplex structure and uses a common set of antennas in the
transmit and receive path. The transmitted carrier and the received
backscatter signal are then separated using circulators. It is important to
note that, unlike the conventional full-duplex wireless links, which transmit
a modulated signal, in BCNs the unmodulated carrier leakage to the receiver
can be efficiently surpassed with low complexity methods [38, 39].

Figure 1: An example of the Reader and BNs structure. The Reader transmits
energy in the downlink, while the BNs transmit data through backscattering the
received energy by changing their variable impedance. During the scheduling
time, the BNs switch to receive mode and the Reader transmits control signals.

The time horizon is slotted to fixed duration time slots. At the beginning of
each time slot, a small portion of the time is devoted to channel estimation
and node discovery. Note that our proposed policy only requires the end-to-end
CSI, $\bm{h}_{n}(t)\bm{h}^{T}_{n}(t)$, which simplifies the CSI estimation
process. The Reader then schedules its transmit and receive beamforming for
the rest of the time slot. Furthermore, the Reader notifies the BNs how much
new data they can admit in their buffers and also sets the values of the
reflection coefficients of the BNs. In the rest of the time slot, the Reader
transmits energy, and concurrently receives the backscattered signal from the
BNs. The channel coefficients are assumed to be constant in a time slot but
vary randomly and independently in consecutive time slots. In time slot $t$,
$h_{n}^{m}(t)\in\mathbb{C},m=1,\ldots,M$, denotes the channel gain of the
reciprocal link between the $m$-th antenna of the Reader and $\text{BN}_{n}$,
and $\bm{h}_{n}(t)=(h_{n}^{1}(t),\ldots,h_{n}^{m}(t))^{T}$ is the channel
coefficients vector associated with $\text{BN}_{n}$.

Figure 2: A time slot structure. At the the beginning of each time slot, a
small portion of time is devoted to node discovery, channel estimation through
pilot transmission and scheduling. In the rest of the time slot the Reader
transmits energy, and the BNs transmit data through backscattering.

###  II-A Downlink Energy Transfer

The Reader continuously transmits a carrier $c(t)$, which is the source of the
energy required for uplink data transmission. Transmission beamforming
technique is adopted in the downlink to increase the energy transmission
efficiency. In time slot $t$, $f_{m}(t)\in\mathbb{C},m=1,\ldots,M$, denotes
the gain of the $m$-th energy transmit path of the Reader, and
$\bm{f}(t)\triangleq\big{(}f_{1}(t),\ldots,f_{M}(t)\big{)}^{T}\in\mathbb{C}^{M}$
denotes the transmit beamforming vector. The carrier $c(t)$ has unit power and
$\lVert\bm{f}(t)\lVert^{2}\leq P$ determines the transmission power of the
Reader, with $P$ denoting the maximum transmission power. Considering the
transmit beamforming vector and the channel gains, the received signal in
$\text{BN}_{n}$ is given by

$\displaystyle e_{n}(t)=c(t)\bm{h}_{n}^{T}(t)\bm{f}(t)+v_{n}(t),$ (1)

where $v_{n}(t)$ is the circular additive white Gaussian noise (AWGN) with
variance $\sigma^{2}_{v}$.

###  II-B Uplink Data Transfer

The BNs modulate the backscattered signal by changing the impedance connected
to the BN’s antenna. In time slot $t$, $x_{n}(t)\in\mathbb{C}$, with
$\operatorname{E}\\{|x_{n}(t)|^{2}\\}=1$, and
$\alpha_{n}(t)\in[0\;\alpha_{\text{max}}]$ denote the modulated signal and the
reflection coefficient of $\text{BN}_{n}$, respectively. Here,
$\alpha_{\text{max}}$ denotes the maximum reflections coefficient of the BNs.
The Reader receives the sum of the backscattered signals from all BNs, that
is,

$\displaystyle z(t)\triangleq\sum_{n\in\mathcal{N}}z_{n}(t)+\bm{w}(t),$ (2)

where, $z_{n}(t)$ is the received signal from $\text{BN}_{n}$, i.e.,

$\displaystyle z_{n}(t)\triangleq\alpha_{n}(t)e_{n}(t)x_{n}(t),$ (3)

and $\bm{w}(t)\in\mathbb{C}^{N}$ denotes the AWGN with covariance matrix
$\sigma_{w}^{2}I_{M}$.

The Reader separates the signal received from different BNs by multiplying the
receive beamforming vector corresponding to each user with the received signal
$z(t)$. Let,
$\bm{g}_{n}(t)\triangleq\big{(}g_{n}^{1}(t),\ldots,g_{n}^{M}(t)\big{)}^{T}\in\mathbb{C}^{M}$,
with $\lVert\bm{g}(t)\rVert^{2}=1$, denote the receive beamforming vector for
$\text{BN}_{n}$, where $g_{n}^{m}(t)\in\mathbb{C},\;\forall m$, is the gain of
the $m$ receive path. The detected signal of $\text{BN}_{n}$ is written as

$\displaystyle
r_{n}(t)\triangleq\bm{g}_{n}^{H}(t)z_{n}(t)+\sum_{\tilde{n}\in\mathcal{N}/n}\bm{g}_{n}^{H}(t)z_{\tilde{n}}(t)+\bm{g}_{n}^{H}(t)\bm{w}(t).$
(4)

Considering the definitions of $z_{n}(t)$ in (3) and $e_{n}(t)$ in (1), the
signal to interference plus noise ratio (SINR) of the signal of
$\text{BN}_{n}$ at the decoder, $\mathcal{S}_{n}(t)$, is given by

$\displaystyle\mathcal{S}_{n}(t)\triangleq\frac{\left|\alpha_{n}(t)\bm{g}_{n}^{H}(t)\bm{h}_{n}(t)\bm{h}_{n}^{T}(t)\bm{f}(t)\right|^{2}}{\sigma^{2}_{w}+\sum_{\tilde{n}\in\mathcal{N}/n}\left|\alpha_{\tilde{n}}(t)\bm{g}_{n}^{H}(t)\bm{h}_{\tilde{n}}(t)\bm{h}_{\tilde{n}}^{T}(t)\bm{f}(t)\right|^{2}}.$
(5)

Note that in practical systems the backscattered noise $v_{n}(t)h_{n}(t)$ is
negligible compared to $w_{n}(t)$ due to channel attenuation [11, 10, 40], and
hence it is neglected in (5). Moreover, the interference due to the
unmodulated carrier leakage from the transmitter can be efficiently surpassed
[14, 9, 11, 10, 38, 39] and, accordingly, it is neglected in the SINR
formulation (5). Considering the SINR of the received backscattered signal and
sufficiently long codewords, the data transmission rate of $\text{BN}_{n}$ in
time slot $t$ is obtained as

$\displaystyle R_{n}(t)=W\log\big{(}1+\mathcal{S}_{n}(t)\big{)},$ (6)

where $W$ is the channel bandwidth (for performance analysis in the cases with
short packets, see Section IV-A).

###  II-C Data Admission

The BNs store the incoming data, e.g., the sensed data from environment, in
their buffers. There is a traffic shaping filter at the input of each buffer
in the BNs which limits the number of input bits in time slot $t$ by
$D_{n}(t)\in[0\;\;D_{\text{max}}],\;\forall n$, with $D_{\text{max}}$ being a
constant. Accordingly, the number of stored data bits in the buffer of
$\text{BN}_{n}$ in time slot $t$, denoted by $Q_{n}(t)$, evolves as

$\displaystyle Q_{n}(t+1)=\big{[}Q_{n}(t)-R_{n}(t)\big{]}^{+}+D_{n}(t).$ (7)

###  II-D Network Controller

There is a network controller at the Reader, which has access to the channel
state information (CSI), $\bm{h}_{n}(t)\bm{h}^{T}_{n}(t)$, and the level of
the stored data in the buffers,
$\bm{Q}(t)\triangleq\big{(}Q_{1}(t),\ldots,Q_{N}(t)\big{)}$. In each time
slot, the controller determines the transmit beamforming vector $\bm{f}(t)$
and the receive beamforming matrix
$\bm{G}(t)\triangleq\big{(}\bm{g}_{1}(t),\ldots,\bm{g}_{N}(t)\big{)}\in\mathbb{C}^{M\times
N}$. Moreover, the controller determines the reflection coefficient vector,
$\bm{\alpha}(t)\triangleq\big{(}\alpha_{1}(t)\ldots,\alpha_{N}(t)\big{)}$, and
the data admission vector,
$\bm{D}(t)\triangleq\big{(}D_{1}(t),\ldots,D_{N}(t)\big{)}$. We design the
network control algorithm to maximize a throughput utility function
$\mathcal{U}(t)$ while keeping the data buffers stable. We consider three
different throughput utility functions, including sum throughput utility

$\displaystyle\mathcal{U}(t)=\mathcal{U}_{s}\big{(}\bm{D}(t)\big{)}\triangleq\sum_{n\in\mathcal{N}}D_{n}(t),$
(8)

proportional throughput utility

$\displaystyle\mathcal{U}(t)=\mathcal{U}_{p}\big{(}\bm{D}(t)\big{)}\triangleq\sum_{n\in\mathcal{N}}\log\big{(}1+D_{n}(t)\big{)},$
(9)

and common throughput utility

$\displaystyle\mathcal{U}(t)=\mathcal{U}_{c}\big{(}\bm{D}(t)\big{)}\triangleq\min_{n\in\mathcal{N}}\big{\\{}D_{n}(t)\big{\\}}.$
(10)

The three introduced utility functions result in different throughput and
fairness tradeoffs. While the sum throughput utility is the most greedy
function in terms of throughput, the common throughput utility is the most
fair utility. The proportional throughput utility, on the other hand, balances
the throughput and fairness. Interestingly, each of these utility functions is
useful in certain applications.

###  II-E Problem Formulation

Let $\mathcal{U}^{\star}$ denote the maximum utility achieved by the optimal
control policy among all policies that stabilize the buffers. In this way,
considering the constraints on the transmission power, reflection coefficient
and data admission, the network controller with maximum utility can be
formulated as the solution of

(11a) (11b) (11c) (11d)

where the expectation in (11) is with respect to the channel randomness.
Constraint (11a) ensures the stability of the buffers. Constraints (11b) and
(11c) are limitations on transmission power and reflection coefficients,
respectively. Moreover, Constraint (11d) determines the maximum number of data
bits that can be admitted to the buffers in each time slot.

Problem (11) is a stochastic utility optimization problem. This problem can be
tackled by dynamic programing (DP) methods. However, DP methods require the
statistical knowledge of the channel state process, which may not be
available. Here, we use the MDPP method [41, Chapter 4] to propose a solution
for Problem (11). MDPP is a general framework for optimizing time averages,
with possibly time average constraints (see, e.g., [42, 43, 44, 45, 46, 47,
48, 49] for different applications of the MDPP method). Using this framework,
a problem with a time average objective function is reduced to a sub-problem
which should be solved in each time slot. Accordingly, following the MDPP
approach, we formulate the optimal $\bm{f}(t),\bm{G}(t),\bm{\alpha}(t)$ and
$\bm{D}(t)$ in each time slot as the solution of an optimization problem with
parameters $\bm{Q}(t)$ and $\bm{h}_{n}(t)\bm{h}_{n}^{T}(t)$. The formulated
problem is non-convex and does not have a trivial solution. However, we use
quadratic and Lagrangian dual transforms [50, 51] to propose an iterative
algorithm for finding the optimal control variables.

##  III The Proposed Control Policy

We follow the MDPP method to propose an online control policy that solves
Problem (11). In summary, we follow these steps:

  1. 1.

We define the Lyapunov function

$\displaystyle L(t)\triangleq\frac{1}{2}\sum_{n\in\mathcal{N}}Q_{n}^{2}(t),$
(11l)

which is a scalar measure of the stored data in all buffers.

  2. 2.

We define the drift-plus-penalty (DPP) function

$\displaystyle\Delta_{p}\big{(}L(t)\big{)}\triangleq\operatorname{E}\Big{\\{}L(t+1)-L(t)\Big{|}\bm{Q}(t)\Big{\\}}-V\operatorname{E}\Big{\\{}\mathcal{U}(t)\Big{\\}},$
(11m)

where $V>0$ is a control parameter. The first term in (11m) shows the drift of
the Lyapunov function in successive time slots. Positive or negative drift
values indicate that the stored data in the buffers have increased or
decreased, respectively. Moreover, the second term is a penalty which
increases as the utility decreases. Intuitively, we expect that under a
control policy, which minimizes (11m) in each time slot, the buffers will be
stable and also the utility will be maximized.

  3. 3.

We introduce an upper bound for the DPP function in Lemma 1. Using the
Lyapunov optimization theorem [41, Theorem 4.2], we show in Theorem
LABEL:th:optimal that under a policy which minimizes the developed upper-
bound, the utility is maximized and the level of the stored data in the
buffers is upper-bounded.

  4. 4.

To find the optimal control policy, we need to minimize the derived upper-
bound in Lemma 1, which includes a non-convex function of
$\bm{f}(t),\bm{G}(t)$ and $\bm{\alpha}(t)$. We use Lagrangian dual and
quadratic transfers [50, 51], and propose an iterative algorithm for finding
the optimal $\bm{f}(t),\bm{G}(t)$ and $\bm{\alpha}(t)$. We also derive the
optimal data admission policy $\bm{D}(t)$.

The details of the analysis are explained as follows. We first introduce an
upper-bound for $\Delta_{p}\big{(}L(t)\big{)}$ in Lemma 1.

######  Lemma 1.

For the DPP function (11m), we have

$\displaystyle\begin{split}\Delta_{p}\big{(}L(t)\big{)}&\leq\Delta_{u}(t)\triangleq
B-\sum_{n\in\mathcal{N}}\operatorname{E}\bigg{\\{}Q_{n}(t)R_{n}(t)\Big{|}\bm{Q}(t)\bigg{\\}}+\\\
&\qquad\qquad\sum_{n\in\mathcal{N}}\operatorname{E}\bigg{\\{}Q_{n}(t)D_{n}(t)\Big{|}\bm{Q}(t)\bigg{\\}}-V\operatorname{E}\bigg{\\{}\mathcal{U}(t)\Big{|}\bm{Q}(t)\bigg{\\}},\end{split}$
(11n)

where

$\displaystyle
B\triangleq\frac{N}{2}\Big{(}D_{\text{max}}+R_{\text{max}}\Big{)}.$ (11o)

Here, $R_{\text{max}}>0$ is a sufficiently large constant, such that we always
have $R_{n}(t)\leq R_{\text{max}}$.

###### Proof.

See Appendix A. ∎

The upper-bound function $\Delta_{u}(t)$ in (11n) is the starting point for
deriving the optimal policy. We formulate the optimization problem

Output: $\bm{f}(t),\bm{g}_{n}(t),\alpha_{n}(t)\;\forall n$.

1:$f(t)\leftarrow\frac{\sum_{n\in\mathcal{N}}Q_{n}(t)\bm{h}^{\ast}_{n}(t)}{\lVert\sum_{n\in\mathcal{N}}Q_{n}(t)\bm{h}^{\ast}_{n}(t)\lVert}\sqrt{P}$.
$\triangleright$ Initialization

2:$\alpha(t)\leftarrow[\alpha_{\text{max}},\ldots,\alpha_{\text{max}}]^{T}$.

3:Calculate $\bm{g}_{n}(t),\;\forall n,$ according to (LABEL:eq:optg).

4:$r_{c}\leftarrow 0$.

5:$\text{it}\leftarrow 0$

6:repeat$\triangleright$ Main loop

7: $r_{p}\leftarrow r_{c}$.

8: Update $\gamma_{n}(t),\;\forall n,$ according to (LABEL:eq:optGamma).

9: Update $y_{n}(t),\;\forall n,$ according to (LABEL:eq:optY).

10: Update $\bm{f}(t)$ according to (LABEL:eq:optf).

11: Update $\alpha_{n}(t),\;\forall n,$ according to (LABEL:eq:optAlpha).

12: Update $\bm{g}_{n}(t),\;\forall n,$ according to (LABEL:eq:optg).

13: Calculate $R_{n}(t),\;\forall n,$ according to (6).

14: $r_{c}\leftarrow\sum_{n\in\mathcal{N}}Q_{n}(t)R_{n}(t)$

15: $\text{it}\leftarrow\text{it}+1$.

16:until $|r_{c}-r_{p}|>\epsilon
r_{p}\land\text{it}\leq\text{it}_{\text{max}}$

###  III-B Data Admission Problem

Considering different utility functions, we derive the corresponding admission
policy through solving Problem (LABEL:prob:DataAdmission).

#### Sum Throughput Utility

Problem (LABEL:prob:DataAdmission) with the sum throughput utility
$\mathcal{U}(t)=U_{s}(\bm{D}(t))$ converts to the linear problem

(11pstxaaaja)

which is solved if

$\displaystyle D_{n}(t)=\begin{cases}D_{\text{max}},&Q_{n}(t)\leq V\\\
0,&\text{Otherwise}.\end{cases}$ (11pstxaaajak)

The policy in (11pstxaaajak) follows a greedy binary rule. Specifically, the
BNs admit the maximum possible data in their buffers, whenever the stored bit
level in the buffer is below $V$, while no data is admitted when the level
reaches $V$.

#### Common Throughput Utility

Adopting the common throughput utility, $U_{c}(\bm{D}(t))$ in Problem
(LABEL:prob:DataAdmission), we have

(11pstxaaajala)

The following Lemma establishes the structure of the solution of Problem
(11pstxaaajal).

######  Lemma 4.

For maximizing the objective function (11pstxaaajal), all BNs should admit
equal amount of data. That is,

$\displaystyle D_{1}(t)=D_{2}(t)=\ldots=D_{N}(t)=D(t).$ (11pstxaaajalam)

###### Proof.

Considering an arbitrary data admission vector $\bm{D}(t)$ that satisfies
(11pstxaaajala), we construct a new vector $\tilde{\bm{D}}(t)$ such that
$\tilde{D}_{1}(t)=\tilde{D}_{2}(t)=\ldots=\min_{n}\\{D_{n}(t)\\}$. Then, we
have

$\displaystyle\sum_{n\in\mathcal{N}}Q_{n}(t)\tilde{D}_{n}(t)-V\min_{n}\\{\tilde{D}_{n}(t)\\}\leq\sum_{n\in\mathcal{N}}Q_{n}(t)D_{n}(t)-V\min_{n}\\{D_{n}(t)\\},$
(11pstxaaajalan)

where equality holds only if $\bm{D}(t)=\tilde{\bm{D}}(t)$. Inequality
(11pstxaaajalan) holds since we have
$\min_{n}\\{D_{n}(t)\\}=\min_{n}\\{\tilde{D}_{n}(t)\\}$ and
$\tilde{D}_{n}(t)\leq D_{n}(t)$. Accordingly, the optimal $\bm{D}(t)$ follows
(11pstxaaajalam). ∎

Using Lemma 4, after adopting $\bm{D}(t)=(D(t),\ldots,D(t))^{T}$, we can
rewrite Problem (11pstxaaajal) as

(11pstxaaajalaoa)

Thus, the optimal data admission rule is

$\displaystyle
D_{n}(t)=\begin{cases}D_{\text{max}},&\sum_{n\in\mathcal{N}}Q_{n}(t)\leq V\\\
0&\text{Otherwise}.\end{cases}$ (11pstxaaajalaoap)

Similar to (11pstxaaajak), the data admission policy in (11pstxaaajalaoap)
follows a binary rule. However, the admission decision in (11pstxaaajalaoap)
is common for all BNs. Accordingly, under the common throughput utility all
BNs will reach the same throughput.

#### Proportional Throughput Utility

Substituting the proportional throughput utility (9) in Problem
(LABEL:prob:DataAdmission), we have

(11pstxaaajalaoaqa)

which is a convex and differentiable function with respect to $\bm{D}(t)$.
Accordingly, comparing the stationary point of (11pstxaaajalaoaq) and the
boundary point in (11pstxaaajalaoaqa), we find the optimal admission policy

$\displaystyle
D_{n}(t)=\begin{cases}D_{\text{max}},&Q_{n}(t)\leq\frac{V}{1+D_{\text{max}}}\\\
0,&Q_{n}(t)\geq V\\\ \frac{V}{Q_{n}(t)}-1,&\text{Otherwise}.\end{cases}$
(11pstxaaajalaoaqar)

The admission policy in (11pstxaaajalaoaqar) is proportional to the inverse of
the level of the stored data in the buffer. This, intuitively, means that the
BN becomes less greedy to admit new data when the buffer size increases.

##  IV Simulation Results

Here, we present the simulation results, and evaluate the performance of our
proposed policy. In all figures, we consider the Rician fading model, that is,

$\displaystyle\bm{h}_{n}(t)=\sqrt{\beta_{n}}\Bigg{(}\sqrt{\frac{K}{K+1}}\bm{h}^{d}_{n}(t)+\sqrt{\frac{1}{K+1}}\bm{h}^{s}_{n}(t)\Bigg{)},$
(11pstxaaajalaoaqas)

where, $\bm{h}^{d}_{n}(t)$ and $\bm{h}^{s}_{n}(t)$ denote the deterministic
and scattered components of the channel, respectively. Moreover, $K$ is the
Rician $K$-factor which determines the ratio between the Rician and the
scattered components, and $\beta_{n}$ represents the path loss factor. Note
that $K=0$ and $K\to\infty$ represent the cases with Rayleigh fading and line-
of-sight channels, respectively. We consider
$\beta_{n}=d_{n}^{-\rho}\left({3\times 10^{8}\over 4\pi f}\right)^{2}$, where
$d_{n}$ is the distance between $\text{BN}_{n}$ and the Reader, $\rho$ is the
path loss exponent and $f$ is the transmit frequency. The entries of the
scattered component vector $\bm{h}^{s}_{n}(t)$ are independent and zero-mean
unit variance circularly symmetric complex Gaussian (CSCG) distributed random
variables. The deterministic components, $\bm{h}^{d}_{n}(t)$, are determined
according to a half wavelength separated uniform linear array setting as
modeled in [55, Eq. 2]. Unless otherwise stated, the BNs are distributed
randomly with a uniform distribution in a circular area around the Reader with
average distance to Reader equal to $30$ m. The general simulation parameters
are as summarized in Table I.

TABLE I: Summary of the simulation parameters. Parameter | Value
---|---
Number of BNs, $N$ | 5
Number of Reader’s antennas, $M$ | 5
Received noise power, $\sigma_{w}$ |  $-110$ dBm
Reader’s maximum transmission power, $P$ |  $500$ mW
Bandwidth, $W$ |  $5$ kHz
Maximum admitted data in each time slot, $D_{\text{max}}$ |  $30$ kbits
Convergence threshold in Algorithm LABEL:alg:ED, $\epsilon$ | 0.01
Maximum iterations in Algorithm LABEL:alg:ED, $\text{it}_{\text{max}}$ | 100
Maximum reflection coefficient, $\alpha_{\text{max}}$ | 0.8
Simulation time, $T$ | $10^{3}$
Rician factor, $K$ | 1
Path loss exponent, $\rho$ | 3
Transmit frequency, $f$ |  $915$ MHz

Considering different parameter settings, Fig. 3 studies the convergence of
Algorithm LABEL:alg:ED. Specifically, in Figures 3(a) and 3(b), the average
number of iterations of Algorithm LABEL:alg:ED is plotted versus the number of
BNs $N$ and the number of Reader’s antennas $M$, respectively. As seen in the
figures, with different parameter settings, Algorithm LABEL:alg:ED converges
with few iterations. However, the average number of iterations increases
slightly as the number of BNs or the number of Reader’s antennas increase.

In Fig. 3(c), we study the effect of the maximum number of iterations
$\text{it}_{\text{max}}$, on the achievable sum throughput utility
$\bar{U}_{s}\triangleq{1\over T}\sum_{t=0}^{T-1}U_{s}(t)$, in a BCN with
$N=\\{5,10\\}$ and $M=\\{5,10\\}$. According to this figure, we achieve the
maximum utility with a few iterations of Algorithm LABEL:alg:ED in each time
slot. Moreover, considering only one iteration, we see in Fig. 3(c) that we
can achieve almost $95\%$ of the maximum utility. Accordingly, we can reduce
the scheduling time through reducing the maximum number of iterations of
Algorithm LABEL:alg:ED significantly, with marginal performance degradation.

(a)

(b)

(c)

Figure 3:  Convergence analysis of Algorithm LABEL:alg:ED with different
parameter settings. Figures (a) and (b) show the average iterations versus $N$
and $M$, respectively. Figure (c) shows the average sum throughput utility
versus the maximum number of iterations $\text{it}_{\text{max}}$.

Considering maximum transmission power $P=\\{100,800\\}$ mW and different
numbers of Reader’s antennas $M=\\{6,8,10\\}$, Fig. 4 demonstrates the
tradeoff between the average proportional throughput utility,
$\bar{U}_{p}={1\over T}\sum_{t=0}^{T-1}U_{p}(t)$, and the average data level
in the buffers. Here, for given values of $P$ and $M$, the utilities are
obtained under different values of $V$ between $10^{7}$ and $10^{9}$. As seen
in Fig. 4, the utility increases as the buffers become more congested.
However, the utility saturates as the average data level in the buffers
increases. This tradeoff between utility and the average data level of the
buffers is in harmony with Theorem LABEL:th:optimal. That is, the optimality
gap is inversely proportional to the maximum data level in the buffers.
Furthermore, considering different parameter settings in Fig. 4, the tradeoff
between utility and the data level in buffers is almost insensitive to the
number of Reader’s antennas $M$ or the transmission power $P$. That is, with
different values of $P$ and $M$ the utility saturates when the data level in
the buffers exceeds $50$ kbits. Finally, as expected, we see in Fig. 4 that
the utility improves as $M$ or $P$ increases. However, this improvement
becomes less dominant as the utility increases.

Figure 4:  Average proportional throughput utility versus the average data
level in the BNs’ buffers for $P=\\{100,800\\}$ mW and $M=\\{6,8,10\\}$.

Considering different values for the average distance of the BNs to the
Reader, Fig. 5 shows the average common throughput utility,
$\bar{U}_{c}={1\over T}\sum_{t=0}^{T-1}U_{c}(t)$, versus the number of
Reader’s antennas $M$. Moreover, the figure compares our scheme with the
proposed policy in [11], which maximizes the minimum achievable rate of all
BNs in each time slot. As can be seen, with the considered average BNs
distances and the number of the Reader’s antennas in Fig. 5, the common
throughput utility under our proposed policy improves on average $13\%$,
compared to the policy in [11]. The reason for such improvement is adopting
buffers in the BNs. The buffers allow the BNs to delay the transmission of the
admitted data in the BNs. Accordingly, the link control policy in Algorithm
LABEL:alg:ED will optimally schedule data transmission for each BN in the best
time slot and, hence, optimizes the resources to maximize the average utility.
Furthermore, Fig. 5 shows the importance of using multi-antenna Readers in
BCNs. According to this figure, considering different average BNs distances,
the utility improves more than $100\%$ when $M$ increases from $6$ to $20$.
Also, adopting more antennas in the Reader, we can increase the coverage
range. For example, consider the parameter settings of Fig. 5 and the common
throughput utility of $1.2$ kbps. Then, by increasing the number of antennas
from $6$ to $10$ the average supported distance increases from 24 m to 36 m.
However, this relative improvement becomes less dominant as the number of
antennas increases.

Figure 5:  Common throughput utility versus the number of the Reader’s
antennas $M$, for the proposed algorithm and the benchmark in [11]. BNs are
distributed randomly in a circular area with average distance to the Reader
$\\{24,30,36\\}$ m and $V=10^{5}$.

In Fig. 6, we compare the average sum throughput utility under our proposed
policy with two other benchmarks. For the first benchmark, called MRT, we
adopt the transmit beamforming according to
$\bm{f}(t)=\sum_{n\in\mathcal{N}}\bm{h}_{n}(t)/\lVert\sum_{n\in\mathcal{N}}\bm{h}_{n}(t)\lVert$
and use the optimal receive beamforming in (LABEL:eq:optg). This specific
choice of $\bm{f}(t)$ is motivated by the fact that MRT beamforming is optimal
for the single BN scenario [54]. For the second benchmark, we consider the
joint transmit and receive beamforming design proposed in [10]. The BNs of the
considered model in [10] have no buffers. Hence, in [10] the sum throughput
utility is optimized in a time slot based framework. Considering
$N\in\\{5,\ldots,12\\}$ BNs distributed uniformly in an area with radius
$\\{50,70\\}$ m, Fig. 6 demonstrates the average sum throughput utility versus
the number of BNs $N$.

As seen in Fig. 6, the sum throughput utility increases almost linearly with
the number of BNs $N$. Moreover, according to this figure, the difference
between the utility under our proposed policy and the policy in [10] is
negligible. Also, both policies provide an average improvement of about $13\%$
and $24\%$ over the MRT with radius equal to $70$ m and $50$ m, respectively.
Accordingly, unlike the common throughput utility, with the sum throughput
utility, adopting buffers in the BNs will not increase the utility. This is
because maximizing the sum throughput utility reduces to opportunistically
maximizing the sum throughput in each time slot. Particularly, considering
$D_{\text{max}}>R_{\text{max}}$, the data admission policy (11pstxaaajak)
implies that during steady state the data buffers in all BNs will fluctuate
near $V+D_{\text{max}}$. Hence, all BNs will have almost equal weights in data
link scheduling problem (LABEL:prob:LinkSchedule) and accordingly, our problem
reduces to maximizing the sum rate as in [10].

Figure 6:  Sum throughput utility versus the number of the BNs $N$ for the
proposed algorithm, the benchmark in [10] and MRT. BNs are distributed
randomly in a circular area with Radius $\\{50,70\\}$ m and $V=10^{7}$.

Considering the utilities, $U_{s}(t)$, $U_{p}(t)$ and $U_{c}(t)$, Fig. 7
compares the average throughput and the received energy of the BNs in a BCN
with 4 BNs. In this BCN, $\text{BN}_{i},i\in\\{1,2,3,4\\}$, are located at
distances 18 m, 22 m, 30 m and 34 m of the Reader, respectively. According to
Fig. 7(a), the sum throughput utility behaves opportunistically where, for
instance, $\text{BN}_{1}$ achieves almost $250\%$ higher throughput compared
to $\text{BN}_{4}$. Whereas, with the common throughput utility all BNs
achieve equal throughputs at a cost of about $20\%$ sum throughput reduction
compared to the other two utilities. The received energy of the BNs in Fig.
7(b) is in harmony with Fig. 7(a). That is, with the sum and proportional
throughput utilities, $\text{BN}_{1}$ receives the highest portion of the
energy transmitted by the Reader. However, with the common throughput utility
the farther BNs receive more energy. Particularly, considering the common
throughput utility, $\text{BN}_{4}$ receives $230\%$ more energy compared to
$\text{BN}_{1}$. While the path loss of the link between $\text{BN}_{4}$ and
the Reader is $(34/18)^{3}\approx 6.7$ times higher than path loss of the link
between $\text{BN}_{1}$ and the Reader.

[width = 1]FairnesComp

(a) Average throughput for each BN.

[width =1]FairnesCompEnergy

(b) Average received energy for each BN.

Figure 7: Utility functions comparison. The distances of the
$\text{BN}_{n},\;n=\\{1,2,3,4\\}$, to the Reader are set to $18$ m, $22$ m,
$30$ m and $34$ m, respectively, and $V=10^{7}$.

###  IV-A Finite Blocklength Analysis of the Proposed Scheme

To simplify the analysis, we presented the results for the cases with
sufficiently long codewords where the Shannon’s capacity, i.e., (6), gives an
appropriate approximation of the achievable rates. However, depending on the
Reader energy budget and the number of nodes, the BCN may be of interest in
the cases with short packets. For this reason, in Fig. 8, we investigate the
effect of the finite length codewords on the communication range under the
proposed policy. Particularly, we define the communication range as the
maximum distance of the BNs to the Reader such that the BNs achieve a
specified common throughput with certain error probability. Here, we consider
BCNs with BNs located circularly around the Reader. Using the fundamental
results of [31, 32, 33, 34] on the achievable rates of short packets, we
replace (6) with [31, Theorem 45],

$\displaystyle R_{n}(t)\approx
W\Bigg{[}\log\left(1+\mathcal{S}_{n}(t)\right)-\sqrt{{1\over{L}}\bigg{[}1-{\bigg{(}1+\mathcal{S}_{n}(t)\bigg{)}^{-2}}\bigg{]}}\mathcal{Q}^{-1}(\psi)\Bigg{]},$
(11pstxaaajalaoaqbh)

which approximates the achievable rate with the finite length codewords. In
(11pstxaaajalaoaqbh), $L$ denotes the packet length, $\mathcal{Q}^{-1}(.)$ is
the inverse $Q$-function and $\psi$ is the maximum codeword error probability
of the decoder. Hence, the second term inside the brackets in
(11pstxaaajalaoaqbh) is notable only in the case of codewords with finite
length. Then, letting $L\rightarrow\infty$ (11pstxaaajalaoaqbh) is simplified
to (6) for the cases with asymptotically long codewords. Note that replacing
(6) with the modified rate function (11pstxaaajalaoaqbh), the results in
Theorem LABEL:th:optimal are still applicable and our data admission policy
remains optimal. However, Algorithm LABEL:alg:ED is an approximate solution
for data link scheduling Problem (LABEL:prob:LinkSchedule) with the modified
rate function (11pstxaaajalaoaqbh).

Considering the common throughputs $\bar{U}_{c}(t)=\\{3,5\\}$ kbps, Fig. 8
shows the communication range versus the numbers of the BNs $N$. This figure
is plotted for codeword lengths $L=\\{\infty,10^{4},10^{3},10^{2}\\}$ and the
error probability $\psi=10^{-3}$. As seen in Fig. 8, in harmony with
intuition, while the Shannon’s capacity-based evaluations are tight for the
cases with long codewords, with finite length codewords the communication
range decreases. This communication range reduction is almost constant for
different number of BNs $N$. Moreover, the communication range reduction is
more significant as the common throughput decreases. This is because, with
lower signal to interference plus noise ratio, the second term inside the
brackets in (11pstxaaajalaoaqbh) is more dominant. Finally, Fig. 8 shows the
tradeoff between communication range and number of BNs. According to this
figure, the communication range decreases almost linearly with the number of
BNs.

Figure 8:  Communication range versus the number of the BNs, $N$ for codeword
lengths $L=\\{\infty,10^{4},10^{3},10^{2}\\}$, target common throughput set to
$\bar{U}_{c}=\\{3,5\\}$ kbps, and $V=10^{7}$.

Considering different codeword lengths $L=\\{10^{2},10^{3}\\}$ and Rician
$K$-factors $K=\\{0,1,10,100\\}$, Fig. 9 shows the average common throughput
utility $\bar{U}_{c}$ versus the maximum error probability $\psi$ in
(11pstxaaajalaoaqbh). According to this figure, considering different
$K$-factors, the utility increases with the error probability almost
logarithmically. However, the common throughput utility is more sensitive to
the error probability with low codeword lengths. As the codeword length
increases, the sensitivity to the error probability decreases. Therefore, at
lower targeted error probabilities short codewords can significantly degrade
the utility which should be considered in BCN design. Moreover, the common
throughput increases as the Rician $K$-factor increases. That is because the
deterministic component of the channel $\bm{h}_{n}^{d}(t)$ becomes more
dominant as $K$ increases. This relative increment saturates at higher
$K$-factors. Finally, the common throughput is more sensitive to the Rician
$K$-factor with longer codewords.

Figure 9:  Common throughput utility versus the maximum codeword error
probability, $\psi$ for codeword lengths $L=\\{10^{2},10^{3}\\}$, Rician
$K$-factors $K\in\\{0,1,10,100\\}$ and $V=10^{7}$.

##  V Conclusion

This paper studied data and energy scheduling in a monostatic BCN with a
multi-antenna Reader and multiple BNs. The BNs adopt buffers to store their
admitted data before transmission to the Reader. We proposed data link control
and data admission policies for maximizing the average value of different
utility functions, including the sum, the proportional and the common
throughput utilities. Through simulation comparisons, we showed the
superiority of our proposed scheme compared to state-of-the-art works.
Specifically, we showed that the sum throughput utility under our proposed
policy achieves the maximum sum channel rate in [10], while stabilizing the
data buffers. Whereas, our proposed policy for optimizing the common
throughput improves the result in [11]. Moreover, using the results on the
achievable rates of finite blocklength codewords, we studied the system
performance in the cases with short packets.

As demonstrated, the proposed policies achieve optimal utility and stabilize
the data buffers in the BNs with few iterations of the data link control
algorithm. Moreover, considering the common throughput utility, adopting
buffers in BNs enables more efficient data scheduling and, hence, improves the
common throughput. Finally, according to our analysis, the finite length of
the codewords affects the communication range significantly, specifically at
low signal to interference plus noise ratios. Also, the utilities are more
sensitive to maximum error probability at short codeword lengths, which should
be carefully compensated for in the BCN design.

##  Appendix A Proof of Lemma 1

Consider the following chain of inequalities,

$\displaystyle\begin{split}\operatorname{E}\Big{\\{}L(t+1)-L(t)\Big{|}\bm{Q}(t)\Big{\\}}&=\frac{1}{2}\sum_{n\in\mathcal{N}}\operatorname{E}\Big{\\{}Q_{n}^{2}(t+1)-Q_{n}^{2}(t)\Big{|}\bm{Q}(t)\Big{\\}}\\\
&\overset{(a)}{=}\frac{1}{2}\sum_{n\in\mathcal{N}}\operatorname{E}\bigg{\\{}\big{(}[Q_{n}(t)-R_{n}(r)]^{+}\big{)}^{2}+D^{2}_{n}(t)+\\\
&\qquad\qquad
2D_{n}(t)[Q_{n}(t)-R_{n}(r)]^{+}-Q_{n}(t)^{2}\Big{|}\bm{Q}(t)\bigg{\\}}\\\
&\overset{(b)}{\leq}\frac{1}{2}\sum_{n\in\mathcal{N}}\operatorname{E}\bigg{\\{}R_{n}(t)^{2}-2Q_{n}(t)R_{n}(t)+D_{n}(t)^{2}+2D_{n}(t)Q_{n}(t)\Big{|}\bm{Q}(t)\bigg{\\}}\\\
&\overset{(c)}{\leq}B+\sum_{n\in\mathcal{N}}\operatorname{E}\bigg{\\{}Q_{n}(t)(D_{n}(t)-R_{n}(t))\Big{|}\bm{Q}(t)\bigg{\\}},\end{split}$
(11pstxaaajalaoaqbi)

where $(a)$ results from (7). Inequality $(b)$ in (11pstxaaajalaoaqbi) holds
because
$\big{(}[Q_{n}(t)-R_{n}(t)]^{+}\big{)}^{2}\leq\big{(}Q_{n}(t)-R_{n}(t)\big{)}^{2}$
and $[Q_{n}(t)-R_{n}(t)]^{+}\leq Q_{n}(t)$. Moreover, inequality $(c)$ comes
from

$\frac{1}{2}\sum_{n\in\mathcal{N}}\big{(}R_{n}^{2}(t)+D_{n}^{2}(t)\big{)}\leq\frac{N}{2}\big{(}R_{max}^{2}+D_{\text{max}}^{2}(t)\big{)}=B.$

Finally, (11n) in Lemma 1 is proved if we subtract
$V\operatorname{E}\big{\\{}\mathcal{U}(t)|\bm{Q}(t)\big{\\}}$ from both sides
of (11pstxaaajalaoaqbi).

##  Appendix B Proof of Theorem LABEL:th:optimal

#### Proof of the first claim

To prove the first claim of the theorem, we show that under the solution of
Problem (11p), we always have $D_{n}(t)=0$ if $Q_{n}(t)\geq V$. To avoid
clutter, we omit the time slot index $t$. Assume that the buffer level in the
$l$-th BN exceeds $V$, i.e., $Q_{l}\geq V$. We define two data admission
vectors $\bm{D}^{1}$ and $\bm{D}^{0}$ that only differ in the $l$-th element.
Specifically, we have $D_{n}^{1}=D_{n}^{0},\;\forall n\neq l$, $D_{l}^{0}=0$
and $D_{l}^{1}>0$. With slightly modified notation, we define
$\Delta_{u}^{i}\triangleq\Delta_{u}\big{(}\bm{f},\bm{G},\bm{\alpha},\bm{D}^{i}\big{)},\;i=\\{0,1\\}$.
Note that we have modified the notation to emphasize that $\Delta_{u}(.)$ in
(11n) is a function $\bm{f},\bm{G},\bm{\alpha}$ and $\bm{D}$. Accordingly, we
have

$\displaystyle\begin{split}\Delta_{u}^{1}-\Delta_{u}^{0}&\overset{(a)}{=}Q_{l}(t)D^{1}_{l}-V\Big{(}U_{\text{g}}\big{(}\bm{D}^{1}\big{)}-U_{\text{g}}\big{(}\bm{D}^{0}\big{)}\Big{)}\\\
&\overset{(b)}{\geq}Q_{l}(t)D^{1}_{l}-VD^{1}_{l}=\big{(}Q_{l}(t)-V\big{)}D^{1}_{l}\geq
0,\end{split}$ (11pstxaaajalaoaqbj)

where, $\text{g}\in\\{\text{s},\text{p},\text{c}\\}$ indicates the sum,
proportional and common throughput utility functions. The equality $(a)$
follows from the definition of $\Delta_{u}(.)$ in (11n). Moreover, $(b)$
results from the fact that the utility functions
$U_{\text{g}}(\bm{D}),\;\text{g}\in\\{\text{s},\text{p},\text{c}\\}$, are
Lipschitz continuous, such that we have

$\displaystyle
U_{\text{g}}\big{(}\bm{D}^{1}\big{)}-U_{\text{g}}\big{(}\bm{D}^{0}\big{)}\leq\sum_{n\in\mathcal{N}}\big{(}D^{1}_{n}-D_{n}^{0}\big{)}=D^{1}_{l}.$
(11pstxaaajalaoaqbk)

Then, (11pstxaaajalaoaqbj) implies that $\Delta_{u}^{1}\geq\Delta_{u}^{0}$,
and hence, a data admission vector with $D_{l}>0$ is not optimal. Accordingly,
we conclude that $Q_{n}(t)$ will never exceed $V+D_{\text{max}}$.

#### Proof of the second claim

The second claim can be proved following the Lyapunov optimization method in
[41, Chapter 4]. According to [41, Theorem 4.5], there is an optimal
stationary policy, denoted by random-only-policy, which is only a function of
$\bm{h}_{n}(t),\;\forall n$. Under the random-only-policy, we have
$\operatorname{E}\Big{\\{}\mathcal{U}(t)|\bm{Q}(t)\Big{\\}}=\mathcal{U}^{\star}$
and $\operatorname{E}\Big{\\{}Q_{n}(t)(D_{n}(t)-R_{n}(t))\Big{\\}}\leq 0$.
Plugging the random-only-policy in $\Delta_{u}(t)$ and using (11n), we have

$\displaystyle\Delta_{p}(L(t))\leq B-V\mathcal{U}^{\star},$
(11pstxaaajalaoaqbl)

where $\Delta_{p}(L(t))$ is evaluated under the solution of Problem (11p).
Note that (11pstxaaajalaoaqbl) holds since the solution of Problem (11p)
minimizes $\Delta_{u}(t)$, and hence, $\Delta_{p}(L(t))$ under the solution of
Problem (11p) is not greater than $\Delta_{u}(t)$ under all alternative
solutions, including the random-only-policy. Averaging both sides of
(11pstxaaajalaoaqbl) over $t=0,\ldots,T-1$, we obtain

$\displaystyle{1\over T}\sum_{n\in\mathcal{N}}Q^{2}_{n}(t)-{1\over
T}\sum_{n\in\mathcal{N}}Q^{2}_{n}(0)-V{1\over
T}\sum_{t=0}^{T-1}\operatorname{E}\\{\mathcal{U}(t)\\}\leq
B-V\mathcal{U}^{\star},$ (11pstxaaajalaoaqbm)

and, rearranging the terms in (11pstxaaajalaoaqbm) and taking
$\lim_{T\rightarrow\infty}$, we have

$\displaystyle\lim_{T\rightarrow\infty}{1\over
T}\sum_{t=0}^{T-1}\operatorname{E}\\{\mathcal{U}(t)\\}\geq\mathcal{U}^{\star}-{B\over
V}.$ (11pstxaaajalaoaqbn)

##  Appendix C Proof of Proposition LABEL:prop:conv

We first show that the value of the function $\sum Q_{n}(t)R_{n}(t)$ is non-
decreasing in each iteration. To avoid clutter, we omit the time slot index
$t$. Moreover, to emphasize the dependence on $\bm{f},\bm{G}$ and
$\bm{\alpha}$, we use the notations
$\tilde{R}(\bm{f},\bm{G},\bm{\alpha},\bm{\gamma})$,
$\tilde{\raisebox{0.0pt}[0.85pt]{$\tilde{R}$}}(\bm{f},\bm{G},\bm{\alpha},\bm{\gamma},\bm{y})$
and $R_{q}(\bm{f},\bm{G},\bm{\alpha})$ instead of $\tilde{R}(t,\bm{\gamma})$,
$\tilde{\raisebox{0.0pt}[0.85pt]{$\tilde{R}$}}(t,\bm{\gamma},\bm{y})$ and
$\sum Q_{n}(t)R_{n}(t)$, respectively. Let $\bm{f}^{i},\bm{G}^{i}$ and
$\bm{\alpha}^{i}$ denote the value of $\bm{f},\bm{G}$ and $\bm{\alpha}$ at the
beginning of the $i$-th iteration. Then, we have

$\displaystyle\begin{split}R_{q}(\bm{f}^{i},\bm{G}^{i},\bm{\alpha}^{i})&\overset{(a)}{=}\tilde{R}(\bm{f}^{i},\bm{G}^{i},\bm{\alpha}^{i},\bm{\gamma}^{\star})\overset{(b)}{=}\tilde{\raisebox{0.0pt}[0.85pt]{$\tilde{R}$}}(\bm{f}^{i},\bm{G}^{i},\bm{\alpha}^{i},\bm{\gamma}^{\star},\bm{y}^{\star})\\\
&\overset{(c)}{\leq}\tilde{\raisebox{0.0pt}[0.85pt]{$\tilde{R}$}}(\bm{f}^{i+1},\bm{G}^{i},\bm{\alpha}^{i},\bm{\gamma}^{\star},\bm{y}^{\star})\overset{(d)}{\leq}\tilde{\raisebox{0.0pt}[0.85pt]{$\tilde{R}$}}(\bm{f}^{i+1},\bm{G}^{i},\bm{\alpha}^{i+1},\bm{\gamma}^{\star},\bm{y}^{\star})\\\
&\leq\max_{\bm{y}}\tilde{\raisebox{0.0pt}[0.85pt]{$\tilde{R}$}}(\bm{f}^{i+1},\bm{G}^{i},\bm{\alpha}^{i+1},\bm{\gamma}^{\star},\bm{y})\overset{(e)}{=}\tilde{R}(\bm{f}^{i+1},\bm{G}^{i},\bm{\alpha}^{i+1},\bm{\gamma}^{\star})\\\
&\leq\max_{\bm{\gamma}}\tilde{R}(\bm{f}^{i+1},\bm{G}^{i},\bm{\alpha}^{i+1},\bm{\gamma})\overset{(f)}{=}R_{q}(\bm{f}^{i+1},\bm{G}^{i},\bm{\alpha}^{i+1})\\\
&\overset{(g)}{\leq}R_{q}(\bm{f}^{i+1},\bm{G}^{i+1},\bm{\alpha}^{i+1}),\end{split}$
(11pstxaaajalaoaqbo)

where $(a)$ holds because $\bm{\gamma}^{\star}$ maximizes
$\tilde{R}(\bm{f}^{i},\bm{G}^{i},\bm{\alpha}^{i},\bm{\gamma})$ and the maximum
value reduces to $R_{q}(\bm{f}^{i},\bm{G}^{i},\bm{\alpha}^{i})$. Similarly,
$(b)$ follows from the fact that $\bm{y}^{\star}$ maximizes
$\tilde{\raisebox{0.0pt}[0.85pt]{$\tilde{R}$}}(\bm{f}^{i},\bm{G}^{i},\bm{\alpha}^{i},\bm{\gamma}^{\star},\bm{y})$,
and the maximum value reduces to
$\tilde{R}(\bm{f}^{i},\bm{G}^{i},\bm{\alpha}^{i},\bm{\gamma}^{\star})$.
Inequalities $(c)$ and $(d)$ hold since $\bm{f}^{i+1}$ and $\bm{\alpha}^{i+1}$
maximize $\tilde{\raisebox{0.0pt}[0.85pt]{$\tilde{R}$}}(.)$ while the other
variables are fixed. The equalities $(e)$ and $(f)$ are concluded similar to
$(b)$ and $(a)$, respectively. Finally, $(g)$ follows since $\bm{G}^{i+1}$
maximizes $R_{q}(\bm{f}^{i+1},\bm{G},\bm{\alpha}^{i+1})$. According to
(11pstxaaajalaoaqbo), the value of the sum $\sum Q_{n}(t)R_{n}(t)$ is non-
decreasing over successive iterations. Moreover, since $\sum Q_{n}(t)R_{n}(t)$
is upper-bounded, the iterations converge.

## References

  * [1] A. Zaidi, A. Bränneby, A. Nazari, M. Hogan, and C. Kuhlins, “Cellular IoT in the 5G era,” Ericsson, White Paper, Feb. 2020.
  * [2] F. Rezaei, C. Tellambura, and S. Herath, “Large-scale wireless-powered networks with backscatter communications–a comprehensive survey,” _IEEE Open J. Commun. Society_ , vol. 1, pp. 1100–1130, Jul. 2020.
  * [3] Y. Zeng, B. Clerckx, and R. Zhang, “Communications and signals design for wireless power transmission,” _IEEE Trans. Commun._ , vol. 65, no. 5, pp. 2264–2290, Mar. 2017.
  * [4] K. Han and K. Huang, “Wirelessly powered backscatter communication networks: Modeling, coverage, and capacity,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 4, pp. 2548–2561, Mar. 2017.
  * [5] W. Liu, K. Huang, X. Zhou, and S. Durrani, “Next generation backscatter communication: systems, techniques, and applications,” _EURASIP J. on Wireless Commun. and Netw._ , vol. 2019, no. 1, p. 69, Mar. 2019.
  * [6] J. Kimionis, A. Bletsas, and J. N. Sahalos, “Increased range bistatic scatter radio,” _IEEE Trans. Commun._ , vol. 62, no. 3, pp. 1091–1104, Feb. 2014.
  * [7] N. Van Huynh, D. T. Hoang, X. Lu, D. Niyato, P. Wang, and D. I. Kim, “Ambient backscatter communications: A contemporary survey,” _IEEE Commun. Surveys Tuts._ , vol. 20, no. 4, pp. 2889–2922, May 2018.
  * [8] G. Yang, C. K. Ho, and Y. L. Guan, “Multi-antenna wireless energy transfer for backscatter communication systems,” _IEEE J. Sel. Areas Commun._ , vol. 33, no. 12, pp. 2974–2987, Sep. 2015.
  * [9] S. Chen, S. Zhong, S. Yang, and X. Wang, “A multiantenna RFID reader with blind adaptive beamforming,” _IEEE Internet Things J._ , vol. 3, no. 6, pp. 986–996, Mar. 2016.
  * [10] D. Mishra and E. G. Larsson, “Sum throughput maximization in multi-tag backscattering to multiantenna reader,” _IEEE Trans. Commun._ , vol. 67, no. 8, pp. 5689–5705, Apr. 2019.
  * [11] ——, “Multi-tag backscattering to mimo reader: Channel estimation and throughput fairness,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 12, pp. 5584–5599, Sep. 2019.
  * [12] G. Yang, Q. Zhang, and Y. Liang, “Cooperative ambient backscatter communications for green internet-of-things,” _IEEE Internet Things J._ , vol. 5, no. 2, pp. 1116–1130, Mar. 2018.
  * [13] S. Idrees, X. Zhou, S. Durrani, and D. Niyato, “Design of ambient backscatter training for wireless power transfer,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 10, pp. 6316–6330, Jun. 2020.
  * [14] D. Li and Y. Liang, “Adaptive ambient backscatter communication systems with MRC,” _IEEE Trans. Veh. Technol._ , vol. 67, no. 12, pp. 12 352–12 357, Sep. 2018.
  * [15] Q. Tao, C. Zhong, X. Chen, H. Lin, and Z. Zhang, “Maximum-eigenvalue detector for multiple antenna ambient backscatter communication systems,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 12, pp. 12 411–12 415, Oct. 2019\.
  * [16] C. Chen, G. Wang, P. D. Diamantoulakis, R. He, G. K. Karagiannidis, and C. Tellambura, “Signal detection and optimal antenna selection for ambient backscatter communications with multi-antenna tags,” _IEEE Trans. Commun._ , vol. 68, no. 1, pp. 466–479, Oct. 2020.
  * [17] C. Chen, G. Wang, H. Guan, Y. Liang, and C. Tellambura, “Transceiver design and signal detection in backscatter communication systems with multiple-antenna tags,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 5, pp. 3273–3288, Feb. 2020.
  * [18] C. Boyer and S. Roy, “Backscatter communication and RFID: Coding, energy, and MIMO analysis,” _IEEE Trans. Commun._ , vol. 62, no. 3, pp. 770–785, Dec. 2014.
  * [19] C. He, S. Chen, H. Luan, X. Chen, and Z. J. Wang, “Monostatic MIMO backscatter communications,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 8, pp. 1896–1909, Jun. 2020.
  * [20] S. Xiao, H. Guo, and Y. Liang, “Resource allocation for full-duplex-enabled cognitive backscatter networks,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 6, pp. 3222–3235, Apr. 2019.
  * [21] Y. Ye, L. Shi, R. Qingyang Hu, and G. Lu, “Energy-efficient resource allocation for wirelessly powered backscatter communications,” _IEEE Commun. Lett._ , vol. 23, no. 8, pp. 1418–1422, Jun. 2019.
  * [22] H. Yang, Y. Ye, and X. Chu, “Max-min energy-efficient resource allocation for wireless powered backscatter networks,” _IEEE Wireless Commun. Lett_ , vol. 9, no. 5, pp. 688–692, Jan. 2020.
  * [23] I. Krikidis, “Retrodirective large antenna energy beamforming in backscatter multi-user networks,” _IEEE Wireless Commun. Lett_ , vol. 7, no. 4, pp. 678–681, Feb. 2018.
  * [24] G. Yang, D. Yuan, Y. Liang, R. Zhang, and V. C. M. Leung, “Optimal resource allocation in full-duplex ambient backscatter communication networks for wireless-powered IoT,” _IEEE Internet Things J._ , vol. 6, no. 2, pp. 2612–2625, Sep. 2019.
  * [25] S. Gong, X. Huang, J. Xu, W. Liu, P. Wang, and D. Niyato, “Backscatter relay communications powered by wireless energy beamforming,” _IEEE Trans. Commun._ , vol. 66, no. 7, pp. 3187–3200, Feb. 2018.
  * [26] X. Wen, S. Bi, X. Lin, L. Yuan, and J. Wang, “Throughput maximization for ambient backscatter communication: A reinforcement learning approach,” in _Proc. IEEE TNEC_ , Chengdu, China, Mar. 2019, pp. 997–1003.
  * [27] N. Van Huynh, D. T. Hoang, D. N. Nguyen, E. Dutkiewicz, D. Niyato, and P. Wang, “Optimal and low-complexity dynamic spectrum access for RF-powered ambient backscatter system with online reinforcement learning,” _IEEE Trans. Commun._ , vol. 67, no. 8, pp. 5736–5752, Apr. 2019.
  * [28] T. T. Anh, N. C. Luong, D. Niyato, Y. Liang, and D. I. Kim, “Deep reinforcement learning for time scheduling in RF-powered backscatter cognitive radio networks,” in _Proc. IEEE WCNC_ , Marrakesh, Morocco, Oct. 2019, pp. 1–7.
  * [29] D. T. Hoang, D. Niyato, P. Wang, D. I. Kim, and L. Bao Le, “Optimal data scheduling and admission control for backscatter sensor networks,” _IEEE Trans. Commun._ , vol. 65, no. 5, pp. 2062–2077, Feb. 2017.
  * [30] T. Liu, X. Qu, and W. Tan, “Online optimal control for wireless cooperative transmission by ambient RF powered sensors,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 9, pp. 6007–6019, Jun. 2020.
  * [31] Y. Polyanskiy, H. V. Poor, and S. Verdu, “Channel coding rate in the finite blocklength regime,” _IEEE Trans. Inf. Theory_ , vol. 56, no. 5, pp. 2307–2359, Jan. 2010.
  * [32] B. Makki, T. Svensson, and M. Zorzi, “Finite block-length analysis of the incremental redundancy HARQ,” _IEEE Wireless Commun. Lett._ , vol. 3, no. 5, pp. 529–532, Aug. 2014.
  * [33] ——, “Wireless energy and information transmission using feedback: Infinite and finite block-length analysis,” _EEE Trans. Commun._ , vol. 64, no. 12, pp. 5304–5318, Sep. 2016.
  * [34] M. Haghifam, M. Robat Mili, B. Makki, M. Nasiri-Kenari, and T. Svensson, “Joint sum rate and error probability optimization: Finite blocklength analysis,” _IEEE Wireless Commun. Lett._ , vol. 6, no. 6, pp. 726–729, Dec. 2017.
  * [35] G. Vannucci, A. Bletsas, and D. Leigh, “A software-defined radio system for backscatter sensor networks,” _IEEE Trans. Wireless Commun._ , vol. 7, no. 6, pp. 2170–2179, Jun. 2008.
  * [36] X. Lu, D. Niyato, H. Jiang, D. I. Kim, Y. Xiao, and Z. Han, “Ambient backscatter assisted wireless powered communications,” _IEEE Wireless Commun._ , vol. 25, no. 2, pp. 170–177, Jan. 2018.
  * [37] N. Fasarakis-Hilliard, P. N. Alevizos, and A. Bletsas, “Coherent detection and channel coding for bistatic scatter radio sensor networking,” _IEEE Trans. Commun._ , vol. 63, no. 5, pp. 1798–1810, Mar. 2015.
  * [38] D. P. Villame and J. S. Marciano, “Carrier suppression locked loop mechanism for UHF RFID readers,” in _Proc. IEEE RFID_ , Orlando, FL, USA, Apr. 2010, pp. 141–145.
  * [39] X. Hao, H. Zhang, Z. Shen, Z. Liu, L. Zhang, H. Jiang, J. Liu, and H. Liao, “A 43.2 $\mu$W 2.4 GHz 64-QAM pseudo-backscatter modulator based on integrated directional coupler,” in _Proc. IEEE ISCAS_ , Florence, Italy, May 2018, pp. 1–5.
  * [40] B. Lyu, D. T. Hoang, and Z. Yang, “User cooperation in wireless-powered backscatter communication networks,” _IEEE Wireless Commun. Lett_ , vol. 8, no. 2, pp. 632–635, Jan. 2019.
  * [41] M. J. Neely, _Stochastic Network Optimization with Application to Communication and Queueing Systems_ , J. Walrand, Ed. Morgan & Claypool Publishers, 2010.
  * [42] L. Huang and M. J. Neely, “Utility optimal scheduling in energy-harvesting networks,” _IEEE/ACM Trans. Netw._ , vol. 21, no. 4, pp. 1117–1130, Dec. 2013.
  * [43] M. J. Neely, “Energy optimal control for time-varying wireless networks,” _IEEE Trans. Inf. Theory_ , vol. 52, no. 7, pp. 2915–2934, Jul. 2006.
  * [44] M. Movahednasab, B. Makki, N. Omidvar, M. R. Pakravan, T. Svensson, and M. Zorzi, “An energy-efficient controller for wirelessly-powered communication networks,” _IEEE Trans. Commun._ , vol. 68, no. 8, pp. 4986–5002, May 2020.
  * [45] M. Movahednasab, N. Omidvar, M. R. Pakravan, and T. Svensson, “Joint data routing and power scheduling for wireless powered communication networks,” in _Proc. IEEE ICC_ , Shanghai, China, May 2019, pp. 1–7.
  * [46] R. Rezaei, N. Omidvar, M. Movahednasab, M. R. Pakravan, S. Sun, and Y. L. Guan, “Efficient, fair, and QoS-aware policies for wirelessly powered communication networks,” _IEEE Trans. Commun._ , vol. 68, no. 9, pp. 5892–5907, Jun. 2020.
  * [47] R. Rezaei, M. Movahednasab, N. Omidvar, and M. R. Pakravan, “Stochastic power control policies for battery-operated wireless power transfer,” in _Proc. IEEE PIMRC_ , Bologna, Italy, Sep. 2018, pp. 1–5.
  * [48] ——, “Optimal and near-optimal policies for wireless power transfer considering fairness,” in _Proc. IEEE GLOBECOM_ , Abu Dhabi, United Arab Emirates, Dec. 2018, pp. 1–7.
  * [49] M. Hadi, M. R. Pakravan, and E. Agrell, “Dynamic resource allocation in metro elastic optical networks using Lyapunov drift optimization,” _J. Opt. Commun. Netw._ , vol. 11, no. 6, pp. 250–259, Jun. 2019.
  * [50] K. Shen and W. Yu, “Fractional programming for communication systems—part I: Power control and beamforming,” _IEEE Trans. Signal Process._ , vol. 66, no. 10, pp. 2616–2630, Mar. 2018.
  * [51] ——, “Fractional programming for communication systems—part II: Uplink scheduling via matching,” _IEEE Trans. Signal Process._ , vol. 66, no. 10, pp. 2631–2644, Mar. 2018.
  * [52] R. W. Freund and F. Jarre, “Solving the sum-of-ratios problem by an interior-point method,” _J. of Global Optimization_ , vol. 19, no. 1, pp. 83–102, Jan. 2001.
  * [53] E. Björnson and E. Jorswieck, _Optimal Resource Allocation in Coordinated Multi-Cell Systems_. Now Foundations and Trends, 2013.
  * [54] E. Björnson, M. Bengtsson, and B. Ottersten, “Optimal multiuser transmit beamforming: A difficult problem with a simple solution structure [lecture notes],” _IEEE Signal Process. Mag._ , vol. 31, no. 4, pp. 142–148, Jun. 2014.
  * [55] Y. Zeng and R. Zhang, “Optimized training design for wireless energy transfer,” _IEEE Trans. Commun._ , vol. 63, no. 2, pp. 536–550, Dec. 2015\.