# Joint Pilot and Payload Power Allocation for Massive-MIMO-enabled URLLC IIoT
Networks

Hong Ren, Cunhua Pan, Yansha Deng, Maged Elkashlan, Arumugam Nallanathan H.
Ren, C. Pan, M. Elkashlan and A. Nallanathan are with School of Electronic
Engineering and Computer Science, Queen Mary University of London, London, E1
4NS, U.K. (Email:h.ren, c.pan, a.nallanathan@qmul.ac.uk). Y. Deng is with the
Department of Informatics, King’s College London, London WC2R 2LS, U.K.
(e-mail:yansha.deng@kcl.ac.uk).

###### Abstract

The Fourth Industrial Revolution (Industrial 4.0) is coming, and this
revolution will fundamentally enhance the way the factories manufacture
products. The conventional wired lines connecting central controller to robots
or actuators will be replaced by wireless communication networks due to its
low cost of maintenance and high deployment flexibility. However, some
critical industrial applications require ultra-high reliability and low
latency communication (URLLC). In this paper, we advocate the adoption of
massive multiple-input multiple output (MIMO) to support the wireless
transmission for industrial applications as it can provide deterministic
communications similar as wired lines thanks to its channel hardening effects.
To reduce the latency, the channel blocklength for packet transmission is
finite, and suffers from transmission rate degradation and decoding error
probability. Thus, conventional resource allocation for massive MIMO
transmission based on Shannon capacity assuming the infinite channel
blocklength is no longer optimal. We first derive the closed-form expression
of lower bound (LB) of achievable uplink data rate for massive MIMO system
with imperfect channel state information (CSI) for both maximum-ratio
combining (MRC) and zero-forcing (ZF) receivers. Then, we propose novel low-
complexity algorithms to solve the achievable data rate maximization problems
by jointly optimizing the pilot and payload transmission power for both MRC
and ZF. Simulation results confirm the rapid convergence speed and performance
advantage over the existing benchmark algorithms.

##  I Introduction

Industry 4.0 has been envisioned as the future paradigm for the next
generation of industrial systems, which integrates advanced manufacturing
functions with the industrial internet-of-things (IIoT) to create a more
intelligent and automatic digital manufacturing system [1]. Traditionally,
industrial control systems mainly rely on wired connections such as cables or
optical fiber, since the current wireless networks cannot meet their stringent
latency and reliability requirements. However, there are some drawbacks to
deploying wired lines. First, significant cost will be incurred by the
installation and maintenance. Second, wired lines are vulnerable to wear and
tear in motion control applications, and suffer from aging. Finally, they
cannot be deployed in some harsh environments, such as those with high
temperatures and rotating part. Hence, to make Industry 4.0 a reality, it is
imperative to design wireless networks tailored for industrial applications to
replace the traditional wired lines. Typical industrial applications require
deterministic communications with ultra reliability ($1-10^{-9}$) and low
latency ($1$ ms), such as factory automation (FA) [2], power system protection
(PSP), and power electronics control (PEC) [3]. Significant research efforts
have been devoted to the design of wireless communications in industrial
applications. However, most of the existing works mainly focused on the
adaption of the upper layers of conventional wireless networks to achieve
deterministic communications, while keeping the physical layer untouched. Some
related standards are WirelessHART, wireless interface for sensors and
actuators (WISA), and Wireless Networks for Industrial Automation/Process
Automation (WIA-PA) [4]. Although keeping the wireless standards of physical
layers can allow faster design and better compatibility, it leads to a
fundamental bottleneck for the system performance. None of the above standards
can meet the stringent demand requested by the most critical FA, PSP and PEC
applications. As a result, more efforts should be devoted to the design from
the physical layer perspective of view.

From the physical layer perspective, the dominate feature of the industrial
applications is that the packet transmission should be completed within short
blocklength due to low latency requirement [5]. Hence, the transmission is not
error-free with any finite/short blocklength channel codes. In this case,
Shannon capacity formula is not applicable since it is based on the principles
of the law of large numbers, and we need to design the resource allocation by
considering the decoding error probability requirement. In [6], Peter _et al._
have derived the approximation formula of the maximum achievable data rate
with finite blocklength transmission, which characterises the complicated
relationships among decoding error probability, channel blocklength, and
signal-to-noise ratio (SNR). Unlike Shannon capacity formula, the data rate
expression under short packet transmission is neither convex nor concave with
respect to the SNR or blocklength [7]. As a result, the optimal resource
allocation under this formula is difficult to obtain.

Recently, there are increasing research studying the transmission design based
on the short packet transmission capacity formula [8, 9, 10, 11, 12, 13, 14].
In specific, the effective throughput maximization was studied in [8] for a
two-device downlink non-orthogonal multiple access (NOMA) system. The overall
error probability is minimized in [9] for a simultaneous wireless information
and power transfer (SWIPT)-enabled decode-and-forward (DF) relaying network.
The decoding error probability minimization was investigated in [10] for a
unmanned aerial vehicle (UAV)-enabled DF relay system. We recently proposed a
low-complexity power and blocklength optimization algorithm for both
orthogonal multiple access (OMA) and NOMA in a two-hop relay system in [11].
However, all these contributions are limited to a simple scenario with two
devices, where one device acts as a relay and the other as the destination
node. In industrial applications, the central controller needs to support a
large number of devices [4]. On one hand, Chen _et al._ in [12] investigated
the effective capacity maximization problem for wireless-powered IoT network
with multiple devices operating in a time division multiple access (TDMA)
mode. However, the time budget is already tight, and the portion allocated to
each device will be marginal. Hence, TDMA strategy is not suitable for the
applications in industrial applications with extremely stringent latency
target. On the other hand, the authors in [13] [14] studied the resource
allocation for multiple devices operating under the orthogonal frequency
division multiple access (OFDMA) mode. Unfortunately, this requires huge
amount of system bandwidth, which is not feasible as some industrial
applications operate over unlicensed spectrum [15].

By equipping a large number of antennas at the base station (BS), massive
multiple-input multiple-output (MIMO) has been widely regarded as the key
enabler for the creation of the fifth generation (5G) wireless networks [16].
By exploiting excessive number of spatial degrees of freedom, massive MIMO is
capable of supporting multiple devices simultaneously without additional time
or frequency resources. In addition, due to the channel hardening effect,
massive MIMO is more immune to the fast fading and can provide deterministic
communications required by the industrial applications. Due to these
attractive advantages, massive MIMO is ideal for supporting industrial
applications with stringent quality of services (QoS) requirements. However,
most of the existing literature adopted Shannon capacity as the performance
metric to optimize the resource allocation [17, 18, 19], which implicitly
assumes the infinite channel blocklength. Therefore, conventional resource
allocation solution based on Shannon capacity is not optimal for industrial
applications with short channel blocklength. To the best of our knowledge, we
are the first to study the resource allocation for massive MIMO providing
ultra-reliability and low-latency communications (URLLC) for any number of
devices. Specifically, our contributions are summarized as follows:

  1. 1.

We derive the closed-form lower bounds (LBs) on the achievable rates for a
uplink massive MIMO system by considering the imperfect channel state
information (CSI) with finite channel blocklength, and assuming both maximum-
ratio combining (MRC) and zero-forcing (ZF) receivers. They can be regarded as
the conventional Shannon capacity minus a penalty term due to short packet
transmission. Simulation results confirm the tightness of the LBs. Given fixed
delay budget, we formulate an optimization problem with the objective to
maximize the weighted sum rate by jointly optimizing the pilot and payload
transmission power subject to the decoding error probability, the minimum data
rates, and the energy constraints for all URLLC devices.

  2. 2.

For the case with MRC receiver, the formulated optimization problem is non-
convex due to the complicated expression of data rate LBs, and it is difficult
to find the globally optimal solution. To deal with this issue, we first
approximate the penalty term in the data LB as a log-function, which can
facilitate the transformation from the original optimization problem to a
series of geometric programs (GPs). Each GP problem can be efficiently solved
with polynomial time. Besides, we provide a novel method to check the
feasibility of the original problem, and provide both complexity and
convergence analysis.

  3. 3.

For the case with ZF receiver, the data rate LB is more complicated and the
algorithm proposed for MRC cannot be directly applied since the nominator of
the signal-to-interference-plus-noise (SINR) is a posynomial function. To
handle this issue, we approximate the polynomial functions with their best
local monomial approximations and then transform the optimization into a
series of GPs with low-complexity. Convergence analysis tailored for the ZF
receiver is further provided.

  4. 4.

Simulation results show that our proposed algorithms converge rapidly, which
verifies the low-complexity of the proposed algorithm. In addition, it is also
shown that the proposed algorithm outperforms the benchmark schemes,
especially the ones adopting Shannon capacity as the optimization performance
metric, which emphasizes the importance of using short packet transmission
theory.

The remainder of this paper is organized as follows. In Section II, system
model and problem formulation are provided. In Section III, we provide a low-
complexity algorithm for joint design of pilot and payload power allocation
for the MRC case. The ZF case is studied in Section IV. Then, simulation
results and analysis are presented in Section V. The final conclusion is drawn
in Section VI.

##  II System Model and Problem Formulation

###  II-A Factory System Model

Consider a uplink multi-device Massive MIMO communication in one factory as
shown in Fig. 1, where the central controller (CC) serves $K$ devices, e.g.,
actuator, robot, etc. The devices need to send their emergency information of
URLLC requirements such as measured data or their current operation states to
the CC. Thus, the CC can process these data information immediately and
provide prompt response/feedback. For simplicity, we focus on the uplink
transmission and the solutions for downlink can be similarly derived. The set
of the devices is denoted as $\mathcal{K}=\\{1,2,\cdots,K\\}$. The CC is
equipped with $M$ antennas and each device is equipped with single antenna due
to their low signal processing capability, where $M\gg K$. Let us denote
$\mathbf{h}_{k}\in\mathbb{C}^{M\times 1}$ as the channel vector from the CC to
the $k$-th device and can be decomposed as
$\textbf{h}_{k}=\sqrt{\alpha_{k}}\bar{\textbf{h}}_{k}$, where $\alpha_{k}$
denotes the large-scale channel gain that includes the pathloss and shadowing,
and $\bar{\textbf{h}}_{k}$ denotes the small-scale fading following the
distribution of $\mathcal{CN}(\textbf{0},\textbf{I})$. Let us denote
$\textbf{H}\in\mathbb{C}^{M\times K}$ as the channel matrix from the $K$
devices to the CC, with
$\textbf{H}=[\textbf{h}_{1},\textbf{h}_{2},\cdots,\textbf{h}_{K}]$.

The $K$ devices need to transmit $K$ packets to the controller. Then, the
$M\times 1$ received signal vector at the CC is given by

$\vspace{-0.15cm}{\bf{y}}=\sum\nolimits_{k\in\mathcal{K}}{{{\bf{h}}_{k}}\sqrt{p_{k}^{d}}}{s_{k}}+{\bf{n}},$
(1)

where $p_{k}^{d}$ is the payload power of the $k$th device, $s_{k}$ is the
zero mean and unit variance Gaussian information message from the $k$th
device, and ${\bf{n}}\sim\mathcal{CN}({\bf{0}},{\bf{I}}_{M})$ is the additive
noise during the data transmission, where the variance of each element is
normalized to unit.

Figure 1: Factory scenario where a massive MIMO central controller serves
multiple devices

###  II-B Channel Estimation in Massive MIMO URLLC

Due to the channel hardening [20] brought by massive MIMO, the system is more
immune to the fast fading, which can provide high reliable services for the
devices. However, to reap the benefits brought by massive MIMO, the CSI should
be available at the CC. Furthermore, TDD mode is always taken as an enabler
for massive MIMO systems since downlink instantaneous CSI is obtained by
estimating uplink CSI based on channel reciprocity [21]. In addition, for
machine type communications, the devices may not be able to perform
complicated signal processing tasks required in frequency division duplexing
(FDD) systems, such as channel estimation calculation, quantization, etc. More
time slots are needed for CSI feedback. Hence, the TDD protocol is adopted in
this paper. All the devices should be allocated with orthogonal pilot
resources so that the CC is able to distinguish the channels from different
devices, thus, the number of symbols for channel estimation should be no
smaller than the number of devices[20].

In a massive MIMO URLLC scenario, each block mainly consists of two parts: 1)
$l_{p}$ symbols for channel estimation (the pilot sequence is of length
$l_{p}$); 2) $l_{d}$ symbols for the $K$ devices’ data transmission, thus the
total number of symbols of the frame is denoted as $L=l_{p}+l_{d}$. As the
symbols used for data transmission is already limited, we assume that $K$
devices sharing the same symbol duration and the frame structure is
illustrated in Fig. 2. Accordingly, the time durations for channel estimation
and data transmission in one frame are given by $t_{p}=l_{p}/B$ and
$t_{d}=l_{d}/B$, respectively, where $B$ is the bandwidth of the system.

Figure 2: Block diagram scheme in TDD URLLC Massive MIMO scheme

In the training phase, all devices simultaneously and synchronously transmit
orthogonal pilot sequences
$\textbf{q}_{1},\cdots,{\textbf{q}}_{l_{p}}\in{\mathbb{C}}^{l_{p}\times 1}$ to
the CC, with ${\textbf{q}_{k}^{H}}\textbf{q}_{k}=1$ and
${\textbf{q}_{i}^{H}}\textbf{q}_{j}=0$, $i\neq j$. Hence, the minimum length
of the pilot sequences to guarantee the orthogonality is equal to $l_{p}=K$.
Based on the received signal, the CC estimates the channel conditions of all
devices, and the received pilot signal at the CC is

$\vspace{-0.2cm}\mathbf{Y}^{p}=\sum\nolimits_{k\in\mathcal{K}}\sqrt{Kp_{k}^{p}}\mathbf{h}_{k}\mathbf{q}_{k}^{H}+\mathbf{N},$
(2)

where $p_{k}^{p}$ is the pilot transmit power at the $k$th device, and
$\mathbf{N}\in\mathbb{C}^{M\times K}$ is the additive Gaussian noise matrix
received during the training phase, whose elements are independently generated
and follow the distribution of $\mathcal{CN}(0,1)$. To obtain channel
$\textbf{h}_{k}$, the CC first multiplies $\textbf{Y}^{p}$ by
$\frac{1}{\sqrt{Kp_{k}^{p}}}\textbf{q}_{k}$, which yields

$\vspace{-0.2cm}\mathbf{y}_{k}^{p}=\frac{1}{\sqrt{Kp_{k}^{p}}}\textbf{Y}^{p}\textbf{q}_{k}=\textbf{h}_{k}+{\textbf{n}}_{k}^{p},$
(3)

where
${\textbf{n}}_{k}^{p}=\frac{1}{\sqrt{Kp_{k}^{p}}}\textbf{N}\textbf{q}_{k}$.
Since $\textbf{q}_{k}$ is a unit-norm vector, it is easy to show that
$\textbf{n}_{k}^{p}$ is still Gaussian distribution whose elements are
independently and identically distributed as
$\mathcal{CN}(0,\frac{1}{Kp_{k}^{p}}{\mathbf{I}}_{M})$. The MMSE estimate of
channel $\textbf{h}_{k}$ is given by

$\hat{\textbf{h}}_{k}=\frac{{{\alpha_{k}}Kp_{k}^{p}}}{{{\alpha_{k}}Kp_{k}^{p}+{1}}}\textbf{y}_{k}^{p},$
(4)

which follows the distribution of
$\mathcal{CN}(\mathbf{0},\sigma_{k}\mathbf{I})$ with $\sigma_{k}$ given by

$\sigma_{k}=\frac{{\alpha_{k}^{2}Kp_{k}^{p}}}{{{\alpha_{k}}Kp_{k}^{p}+1}}.$
(5)

According to the property of MMSE estimation, channel estimation error
$\tilde{\textbf{h}}_{k}=\textbf{h}_{k}-\hat{\textbf{h}}_{k}$ is independent of
$\hat{\textbf{h}}_{k}$, and follows the distribution of
$\mathcal{CN}(\mathbf{0},\delta_{k}\mathbf{I}_{M})$, where $\delta_{k}$ is
given by

$\delta_{k}=\frac{{{\alpha_{k}}}}{{{\alpha_{k}}Kp_{k}^{p}+1}}.$ (6)

###  II-C Achievable Data Rate for Massive MIMO URRLC

By taking into account the number of symbols for pilot transmission, to
achieve the decoding error probability of $\varepsilon_{k}$ for the $k$th
device, the instantaneous achievable data rate $R_{k}$ can be accurately
approximated by [22]

$\vspace{-0.05cm}{R_{k}}\approx{\left({1-\beta}\right){{\log}_{2}}(1+{\gamma_{k}})-\sqrt{\frac{{\left({1-\beta}\right){V_{k}}}}{L}}\frac{{{Q^{-1}}({\varepsilon_{k}})}}{{\ln
2}}},$ (7)

where $\beta$ is equal to
$\beta={{{K}}\mathord{\left/{\vphantom{{{K}}L}}\right.\kern-1.2pt}L}$,
$\gamma_{k}$ is the signal to interference plus noise ratio (SINR) of the
$k$th device, $Q^{-1}$ is the inverse function
$Q(x)=\frac{1}{\sqrt{2\pi}}\int_{x}^{\infty}e^{-\frac{t^{2}}{2}}\text{d}t$,
and $V_{k}$ is the channel dispersion given by $V_{k}=1-(1+\gamma_{k})^{-2}$.
As seen from (7), when the blocklength $L$ approaches infinity, the data rate
$R_{k}$ will approach $\left({1-\beta}\right){\log_{2}}(1+\gamma_{k})$, which
is the classic Shannon capacity. The second term in (7) can be interpreted as
a penalty on the rate in order to guarantee the decoding error probability
$\varepsilon_{k}$.

In the following, we derive the expression of $\gamma_{k}$ for two different
low-complexity detection schemes: 1) maximum-ratio combining (MRC); 2) zero-
forcing (ZF).

Define the estimated channels as
$\widehat{\textbf{H}}=[\hat{\textbf{h}}_{1},\hat{\textbf{h}}_{2},\cdots,\hat{\textbf{h}}_{K}]$
and channel estimation errors as ${\widetilde{\bf H}}=[{{{\tilde{\bf
h}}}_{1}},{{{\tilde{\bf h}}}_{2}},\cdots,{{{\tilde{\bf h}}}_{K}}]$. Let
$\bf{A}$ be an $M\times K$ linear detection matrix that is based on the
estimated channel $\widehat{\textbf{H}}$. By using the linear detection
$\bf{A}$, the received signal can be processed as

$\vspace{-0.3cm}{\bf{y}}^{D}={{\bf{A}}^{H}}{\bf{y}}.$ (8)

Two conventional low-complexity linear detectors are considered:

$\vspace{-0.05cm}{\bf{A}}=\left\\{\begin{array}[]{l}\widehat{\bf{H}},\qquad\qquad\quad{\rm{for}}\quad{\rm{MRC}}\\\
\widehat{\bf{H}}{\left({{{\widehat{\bf{H}}}^{H}}\widehat{\bf{H}}}\right)^{-1}},\
{\rm{for}}\quad{\rm{ZF}}.\end{array}\right.$ (9)

Then, the processed signal after using the detector is given by

$\displaystyle\vspace{-0.35cm}{\bf{y}}^{D}$ $\displaystyle=$
$\displaystyle{{\bf{A}}^{H}}\sum\nolimits_{i\in\mathcal{K}}{{{\bf{h}}_{i}}\sqrt{p_{i}^{d}}}{s_{i}}+{{\bf{A}}^{H}}{\bf{n}}.$
(10) $\displaystyle=$
$\displaystyle{{\bf{A}}^{H}}\sum\nolimits_{i\in\mathcal{K}}{{{{\hat{\bf
h}}}_{i}}\sqrt{p_{i}^{d}}}{s_{i}}+{{\bf{A}}^{H}}\sum\nolimits_{i\in\mathcal{K}}{{{{\tilde{\mathbf{h}}}}_{i}}\sqrt{p_{i}^{d}}}{s_{i}}+{{\bf{A}}^{H}}{\bf{n}},\vspace{-0.2cm}$

where the last equality is obtained by using ${{\bf{h}}_{k}}={{{\hat{\bf
h}}}_{k}}+{{{\tilde{\bf h}}}_{k}}$. The detection signal for the $k$th device
is given by

$\vspace{-0.1cm}{{\bf{y}}^{D}_{k}}={\mathbf{a}}_{k}^{H}{{{\hat{\bf
h}}}_{k}}\sqrt{p_{k}^{d}}{s_{k}}+{\bf{a}}_{k}^{H}\sum\nolimits_{i\in\mathcal{K}\setminus
k}{{{{\hat{\bf
h}}}_{i}}\sqrt{p_{i}^{d}}}{s_{i}}+{\bf{a}}_{k}^{H}\sum\nolimits_{i\in\mathcal{K}}{{{{\tilde{\mathbf{h}}}}_{i}}\sqrt{p_{i}^{d}}}{s_{i}}+{\mathbf{a}}_{k}^{H}{\mathbf{n}}.$
(11)

where ${\bf{a}}_{k}$ is the $k$th column of matrix ${\bf{A}}$. Since
$\widehat{\textbf{H}}$ and $\widetilde{\textbf{H}}$ are independent,
${\bf{a}}_{k}$ is also independent of $\widetilde{\textbf{H}}$. The CC will
treat the estimated channel as the true channel, and the last three terms of
(11) are regarded as interference and noise. Then, the SINR for the $k$th
device $\gamma_{k}$ is given by

${\gamma_{k}}=\frac{{p_{k}^{d}{{\left|{{\bf{a}}_{k}^{H}{{{\hat{\bf
h}}}_{k}}}\right|}^{2}}}}{{\sum\nolimits_{i\in\mathcal{K}\setminus
k}{p_{i}^{d}{{\left|{{\bf{a}}_{k}^{H}{{{\hat{\bf
h}}}_{i}}}\right|}^{2}}}+\sum\nolimits_{i\in\mathcal{K}}{p_{i}^{d}{{\left|{{\bf{a}}_{k}^{H}{{{\tilde{\bf
h}}}_{i}}}\right|}^{2}}}+{{\left\|{{{\bf{a}}_{k}}}\right\|}^{2}}}}.$ (12)

Remark: Massive MIMO can offer the channel hardening effect, where the channel
variations between different channel fading blocks can be averaged out and the
achievable data rate for these fading blocks mainly depend on the large-scale
fading, which changes very slowly. As a result, we first derive the lower
bound (LB) for the achievable data rate as a function of large-scale channel
fading parameters, and then optimize the resource allocation based on the
large-scale fading information rather than the small-scale fading, which can
significantly _reduce the computational delay_ that is beneficial for URLLC
applications. In other words, when the large-scale fading parameters of all
devices are given, we can use our developed algorithm to find the optimal
power allocation, which can be used for consecutive channel fading blocks. The
algorithms are needed to be rerun only when the large-scale channel fading
parameters have changed, which vary much slowly compared with the small-scale
channel fading. To make it more clear, a block diagram scheme is given in Fig.
3, where the large-scale channel gains of any devices vary at $t=t_{0}^{Ls}$
and $t=t_{1}^{Ls}$. In general, the channel coherence time is much longer than
the packet transmission time as we consider the URLLC services, and the time
difference for large-scale fading ($t_{1}^{Ls}-t_{0}^{Ls}$) is much larger
than the channel coherence time. The power allocation obtained at time
$t=t_{0}^{Ls}$ can be employed for the subsequent transmissions until
$t=t_{1}^{Ls}$, which significantly reduce the computational time. Please note
in our scheme, the channel estimation should be performed at the beginning of
each channel coherence since the decoding requires the channel state
information as shown in (9), while the power allocation needs to be updated
once the large-scale fading gains change.

Figure 3: Block diagram for transmission in TDD URLLC Massive MIMO scheme.

Due to the channel hardening effect, in this paper we focus on the ergodic
achievable data rate that is defined as
${\bar{R}_{k}}={\mathbb{E}}\left\\{{R_{k}}\right\\}$, where the expectation is
taken over the randomness of $\left\\{{{{{\hat{\bf h}}}_{k}},{{{\tilde{\bf
h}}}_{k}},\forall k}\right\\}$. Unfortunately, the exact average achievable
data rate ${\bar{R}_{k}}$ with channel uncertainty is not available. In the
following, we aim to derive the closed-form expression of the LB of the rate
expression, which is more tractable to analyse and optimize. To this end, we
first define function $f(x)$ as

$f(x)=\ln\left({1+\frac{1}{x}}\right)-a\sqrt{\frac{{2x+1}}{{{{(x+1)}^{2}}}}}\geq
0,x>0$ (13)

where $a$ is a fixed positive value. In the following, we derive the
feasible region of function $f(x)$. Since $f(x)\geq 0$, from (13) we have

$a\leq\frac{{(x+1)\ln\left({1+\frac{1}{x}}\right)}}{{\sqrt{2x+1}}}\buildrel\Delta\over{=}g(x).$
(14)

The first-order derivative of $g(x)$ with respect to $x$ is given by

$g^{\prime}(x)=\frac{{-2-\frac{1}{x}+x\ln\left({1+\frac{1}{x}}\right)}}{{{{\left({2x+1}\right)}^{\frac{3}{2}}}}}\leq\frac{{-1-\frac{1}{x}}}{{{{\left({2x+1}\right)}^{\frac{3}{2}}}}}<0$
(15)

where the first inequality follows by using the relation
$\ln\left({1+\frac{1}{x}}\right)<\frac{1}{x}$. Hence, $g(x)$ is a
monotonically decreasing function of $x$. In addition,
$\mathop{\lim}\limits_{x\to 0}g(x)=\infty$ and
$\mathop{\lim}\limits_{x\to\infty}g(x)=0$, where the latter equation is
obtained by using the L’Hospital’s rule. Hence, from (14), we know that the
feasible region of $f(x)$ is given by

$\vspace{-0.45cm}{\bf\Omega}=\left\\{{\left.x\right|0<x\leq{g^{-1}}(a)}\right\\}.$
(16)

Then, we have the following lemma.

_Lemma 1_ : Function $f(x)$ defined in (13) is decreasing and convex for
$x\in\bf\Omega$.

_Proof_ : Please refer to Appendix B in our recent work [23].

Based on Lemma 1, we are able to derive the LB of ${\bar{R}_{k}}$ in the
following. The instantaneous data rate $R_{k}$ in (7) can be written as
follows:

$\vspace{-0.05cm}{R_{k}}=\frac{{1-\beta}}{{\ln
2}}{f_{k}}\left({\frac{1}{{{\gamma_{k}}}}}\right),$ (17)

where function ${f_{k}}(\cdot)$ is in the same format as $f(\cdot)$ in (13),
where the parameter $a$ is

$\vspace{-0.3cm}{a_{k}}={{{Q^{-1}}({\varepsilon_{k}})}}/{\sqrt{L\left({1-\beta}\right)}}.$
(18)

In addition, we need to guarantee that $R_{k}\geq 0$, and thus
$\gamma_{k}\geq{1\mathord{\left/{\vphantom{1{{g^{-1}}({a_{k}})}}}\right.\kern-1.2pt}{{g^{-1}}({a_{k}})}}$.
By using lemma 1 and Jensen’s inequality, we obtain the following LB on the
ergodic data rate:

$\vspace{-0.05cm}{{\bar{R}}_{k}}\geq{{\hat{R}}_{k}}\buildrel\Delta\over{=}\frac{{1-\beta}}{{\ln
2}}{f_{k}}\left({{\mathbb{E}}\left\\{{\frac{1}{{{\gamma_{k}}}}}\right\\}}\right).$
(19)

In the following theorem, we derive the expression of ${{\hat{R}}_{k}}$ for
each beamforming solution.

_Theorem 1_ : The ergodic achievable rate for the $k$th device for MRC in
finite blocklenghth regime can be lower bounded by:

$\vspace{-0.1cm}{{\hat{R}}_{k}}\buildrel\Delta\over{=}\frac{{1-\beta}}{{\ln
2}}{f_{k}}\left(\frac{1}{{{{\hat{\gamma}}_{k}}}}\right)$ (20)

where ${\hat{\gamma}}_{k}$ is given by

$\vspace{-0.1cm}{{\hat{\gamma}}_{k}}=\frac{{p_{k}^{d}(M-1){\sigma_{k}}}}{{\sum\nolimits_{i\in\mathcal{K}\setminus
k}{p_{i}^{d}{\sigma_{i}}}+\sum\nolimits_{i\in\mathcal{K}}{p_{i}^{d}{\delta_{i}}}+1}}.$
(21)

_Proof_ : Please refer to Appendix A.

For the ZF detection, the LB of the ergodic data rate is given by the
following theorem.

_Theorem 2_ : The ergodic achievable rate for the $k$th device for ZF in
finite blocklength regime is lower bounded by:

$\vspace{-0.3cm}{{\hat{R}}_{k}}\buildrel\Delta\over{=}\frac{{1-\beta}}{{\ln
2}}{f_{k}}\left(\frac{1}{{{{\hat{\gamma}}_{k}}}}\right)$ (22)

where ${\hat{\gamma}}_{k}$ is given by

$\vspace{-0.2cm}{{\hat{\gamma}}_{k}}=\frac{{(M-K){\sigma_{k}}p_{k}^{d}}}{{\sum\nolimits_{i\in\mathcal{K}}{p_{i}^{d}{\delta_{i}}+1}}}.$
(23)

_Proof_ : Please refer to Appendix B.

In the proof of Theorem 1 and Theorem 2, we notice that the above derived LB
of the ergodic data rates are valid for any number of antennas ($M>K$ for the
ZF). However, the gap between the LB and the actual ergodic data rate is
reduced when the number of antennas is large, which is the case for massive
MIMO system. These LBs are commonly used in the literature concerning massive
MIMO systems. Hence, we use these lower bounds throughout this paper. In the
following, we aim to optimize the pilot power $p_{k}^{p}$ and payload power
$p_{k}^{d}$ to maximize the weighted sum rate. The optimization is performed
only when any large-scale fading parameter changes. The simulation results in
Section V also verify the tightness of the derived LBs. Hence, the
optimization solutions are applicable for _a large time scale_ , which is
appealing for URLLC applications.

###  II-D Problem Formulation

In this paper, we jointly optimize the power allocation for pilot and data
transmission of each device for maximizing the weighted sum rate of all
devices. By using Theorem 1, the weighted sum rate maximization problem can be
formulated as

$\displaystyle\mathop{\max}\limits_{\\{p_{k}^{p}\\},\\{p_{k}^{d}\\}}\;\;\;$
$\displaystyle\sum\nolimits_{k\in\mathcal{K}}{{w_{k}}{{\hat{R}}_{k}}}$ (24a)
$\displaystyle{\rm{s.t.}}\;\;\;$ $\displaystyle{{\hat{R}}_{k}}\geq
R_{k}^{{\rm{req}}},\forall k,$ (24b)
$\displaystyle{K}p_{k}^{p}+(L-K)p_{k}^{d}\leq{E_{k}},\forall k,$ (24c)

where ${{\hat{R}}_{k}}$ and ${{\hat{\gamma}}_{k}}$ are given in (20) and (21)
for MRC and (22) and (23) for ZF, respectively, $w_{k}$ is the weight of
device $k$ used to guarantee the fairness among the devices, constraint (24b)
denotes the minimum data rate requirement for the $k$-th device, constraint
(24c) means the energy constraint for each device.

Power control for weighted sum rate problem with interference is well known to
be an NP-hard problem even under perfect CSI [24]. It becomes more complicated
for the more general case with imperfect CSI and finite blocklength. In this
paper, we aim for designing efficient algorithms with polynomial-time
complexity to solve the weighted sum rate problem with imperfect CSI.

To this end, we first simplify the problem formulation in (24). The first-
order derivative of ${{\hat{R}}_{k}}$ w.r.t. ${\hat{\gamma}}_{k}$ is given by
${{\hat{R}^{\prime}}_{k}}=-\frac{{\left({1-\beta}\right)}}{{\hat{\gamma}_{k}^{2}\ln
2}}{f_{k}^{\prime}}\left({\frac{1}{{{{\hat{\gamma}}_{k}}}}}\right)\geq 0$
111Since ${{\hat{R}}_{k}}\geq R_{k}^{{\rm{req}}}>0$, the feasible
${{\hat{\gamma}}_{k}}$ must lie in the range of
${\bf\Omega}=\left\\{{\left.{\hat{\gamma}}_{k}\right|0<1/{\hat{\gamma}}_{k}\leq{g^{-1}}(a_{k})}\right\\}$.
Hence, Lemma 1 holds and we have
${f_{k}^{\prime}}\left({\frac{1}{{{{\hat{\gamma}}_{k}}}}}\right)\geq 0$..
Hence, constraint (24b) can be transformed as

$\vspace{-0.1cm}{{\hat{\gamma}}_{k}}\geq{1\mathord{\left/{\vphantom{1{{f_{k}^{-1}}\left({\frac{{R_{k}^{{\rm{req}}}\ln
2}}{{1-\beta}}}\right)}}}\right.\kern-1.2pt}{{f_{k}^{-1}}\left({\frac{{R_{k}^{{\rm{req}}}\ln
2}}{{1-\beta}}}\right)}},\forall k.$ (25)

To additionally simplify the problem formulation in (24), we introduce
auxiliary variables ${\chi_{k}},\forall k$, and then Problem (24) can be
equivalently transformed as follows

$\displaystyle\mathop{\max}\limits_{\\{\chi_{k}\\},\\{p_{k}^{p}\\},\\{p_{k}^{d}\\}}\;\;\;$
$\displaystyle\sum\nolimits_{k\in\mathcal{K}}{{\tilde{w}}_{k}}\left[{\ln(1+{\chi_{k}})-{a_{k}}G({\chi_{k}})}\right]$
(26a) $\displaystyle{\rm{s.t.}}\;\;\;$
$\displaystyle{{\hat{\gamma}}_{k}}\geq\chi_{k},\forall k,$ (26b)
$\displaystyle(\ref{fhhthtyh}),(\ref{freocdsdi}),$ (26c)

where ${{\tilde{w}}_{k}}=\frac{{\left({1-\beta}\right){w_{k}}}}{{\ln 2}}$ and
$G(\chi_{k})=\sqrt{1-\frac{1}{{{{(1+\chi_{k})}^{2}}}}}$. The equivalence
between Problem (26) and Problem (24) can be readily proved by using the
contradiction method. Problem (26) and Problem (24) are equivalent in the
sense that they have the same power allocation solutions and the same
objective function (OF) value. Hence, in the following, we focus on solving
Problem (26).

Due to different SINR expressions, in the following two sections, we optimize
the power allocation for MRC and ZF, respectively.

##  III Weighted Sum Data Rate for MRC

In this section, we aim to deal with weighted sum rate maximization problem
for MRC.

###  III-A Algorithm Design

The complicated function $G(\chi_{k})$ in (26a) makes the optimization problem
difficult to solve. To resolve this issue, we first study several properties
of this function.

_Lemma 2_ : $G({x})$ is a concave function of $x$.

_Proof_ : Please refer to Appendix C.

For the URLLC applications, the decoding error probability ${\varepsilon_{k}}$
for each device is much smaller than $0.5$ so that $a_{k}$ is a strictly
positive value, and the OF of Problem (26) is to maximize the difference of
two concave functions. However, due to the non-convex constraints in Problem
(26). This problem does not belong to the class of difference of convex (DC)
problem. In addition, due to the additional term of ${a_{k}}G({\chi_{k}})$ in
the OF of Problem (26), this problem cannot be solved by dealing with a
sequence of geometric programs (GPs) as in [25], which considered Shannon
capacity formula under the infinite blocklength regime. The intuitive method
to solve Problem (26) is to approximate function $G(\chi_{k})$ as its first-
order Taylor series expansion and solve the approximated problem until
convergence. However, the approximation is in linear form while the first term
in the OF is in log-function form. The successive GP method in [25] is still
not applicable. To deal with this issue, we approximate function $G(x)$ as a
log-function as shown in the following lemma.

_Lemma 3_ : For any given $\tilde{x}\geq\frac{{\sqrt{17}-3}}{4}$, the
following inequality holds:

$G(x)=\sqrt{1-\frac{1}{{{{(1+x)}^{2}}}}}\leq\rho{\rm{ln}}(x)+\eta\buildrel\Delta\over{=}F(x),\forall
x\geq\frac{{\sqrt{17}-3}}{4},$ (27)

where $\rho$ and $\eta$ are given by

${\rho}=\frac{{\tilde{x}}}{{\sqrt{{{\tilde{x}}^{2}}+2\tilde{x}}}}-\frac{{\tilde{x}\sqrt{{{\tilde{x}}^{2}}+2\tilde{x}}}}{{{{\left({1+\tilde{x}}\right)}^{2}}}},$
(28)

and

$\vspace{-0.2cm}{\eta}=\sqrt{1-\frac{1}{{{{\left({1+\tilde{x}}\right)}^{2}}}}}-{\rho}\ln(\tilde{x}).$
(29)

In addition, we have:

$\vspace{-0.5cm}G(\tilde{x})=F(\tilde{x}),G^{\prime}(\tilde{x})=F^{\prime}(\tilde{x}),$
(30)

which means that the approximation $F(x)$ is tight at $x=\tilde{x}$.

_Proof_ : Please refer to Appendix D.

According to (16), the optimization variable ${\hat{\gamma}}_{k}$ should be no
smaller than
${1\mathord{\left/{\vphantom{1{{g^{-1}}({a_{k}})}}}\right.\kern-1.2pt}{{g^{-1}}({a_{k}})}}$.
Note that
${1\mathord{\left/{\vphantom{1{{g^{-1}}({a_{k}})}}}\right.\kern-1.2pt}{{g^{-1}}({a_{k}})}}$
is a decreasing function of the decoding error probability and blocklength,
while an increasing function of the number of device. Hence, for typical FA
cases [26], where the typical required error probability is lower than
${10^{-8}}$, the available channel blocklength is smaller than 200, and the
number of devices that should be supported is larger than 5,
${1\mathord{\left/{\vphantom{1{{g^{-1}}({a_{k}})}}}\right.\kern-1.2pt}{{g^{-1}}({a_{k}})}}$
is larger than $\frac{{\sqrt{17}-3}}{4}$, which implies that Lemma 3 is
applicable for our considered optimization problem in (26).

In Fig. 2, we compare the approximation accuracy of function $F(x)$ and the
linear approximation function $S(x)$ defined as follows:

$\vspace{-0.1cm}G(x)\leq
G(\tilde{x})+G^{\prime}(\tilde{x})\left({x-\tilde{x}}\right)\buildrel\Delta\over{=}S(x),$
(31)

where the inequality holds since $G(x)$ is a concave function proved in
Lemma 2. It can be observed from Fig. 4 that the approximation function $F(x)$
is more accurate than the linear function $S(x)$ over the whole region of $x$
at different points of $\tilde{x}$. We can also find that the curve of the
log-function $F(x)$ is always above the curve of the linear function $S(x)$,
which verifies the correctness of the theoretical conclusion in Lemma 3. The
most important advantage of approximating $G(x)$ as the log-function $F(x)$ is
that we can transform the original problem into a GP problems that enables us
to find the optimal solution.

Figure 4: Approximation error for different approximation functions at three
different points: a) $\tilde{x}=0.5$; b) $\tilde{x}=3$; c) $\tilde{x}=6$.

In the following lemma, we also provide the LB of $\ln(1+{\chi_{k}})$, which
enables us to develop low-complexity algorithms.

_Lemma 4_ : For any given $\tilde{x}\geq 0$, function $\ln(1+x)$ is lower
bounded by

$\vspace{-0.3cm}\ln(1+x)\geq{\hat{\rho}}\ln x+\hat{\eta},\forall x\geq 0,$
(32)

where $\hat{\rho}$ and $\hat{\eta}$ are given by

$\vspace{-0.3cm}\hat{\rho}=\frac{{\tilde{x}}}{{1+\tilde{x}}},\hat{\eta}=\ln(1+\tilde{x})-\frac{{\tilde{x}}}{{1+\tilde{x}}}\ln\tilde{x}.$
(33)

In addition, the bound is tight at $x=\tilde{x}$.

_Proof_ : The proof is similar to those in Lemma 3, and thus omitted for
simplicity.

Based on Lemma 3 and Lemma 4, we are now ready to solve Problem (26). The main
idea is to approximate the OF of Problem (26) as the approximated functions
provided in Lemma 3 and Lemma 4, and then solve the approximate problem in an
iterative manner. In the following, we provide the details of the iterative
algorithm.

First, we denote ${\bf{P}}^{(i)}=\\{p_{k}^{p(i)},p_{k}^{d(i)},\forall k\\}$ as
the power allocation in the $i$-th iteration, and the corresponding $\chi_{k}$
is given by $\chi_{k}^{(i)}$. Then, in the $i+1$-th iteration, we can
approximate $G(\chi_{k})$ around $\chi_{k}^{(i)}$ as function
$F(\chi_{k})=\rho_{k}^{(i)}{\rm{ln}}\chi_{k}+\eta_{k}^{(i)}$, where
$\rho_{k}^{(i)}$ and $\eta_{k}^{(i)}$ are obtained from (28) and (29)
respectively with $\tilde{x}=\chi_{k}^{(i)}$. In addition, we can approximate
$\ln(1+\chi_{k})$ around $\chi_{k}^{(i)}$ as
${\hat{\rho}_{k}^{(i)}}\ln\chi_{k}+{\hat{\eta}_{k}^{(i)}}$, where
${\hat{\rho}_{k}^{(i)}}$ and ${\hat{\eta}_{k}^{(i)}}$ can be obtained from
(33) with $\tilde{x}=\chi_{k}^{(i)}$. By substituting these approximations
into (26a) and recalling that $a_{k}$ is a positive value, we can obtain the
LB of the OF by using Lemma 3 and Lemma 4 as follows

$\sum\nolimits_{k\in\mathcal{K}}{{\tilde{w}}_{k}}\left[{\ln(1+{\chi_{k}})-{a_{k}}G({\chi_{k}})}\right]\geq\sum\nolimits_{\mathcal{K}}{{\tilde{w}}_{k}}\left[{{\hat{\rho}_{k}^{(i)}}\ln\chi_{k}+{\hat{\eta}_{k}^{(i)}}-a_{k}\rho_{k}^{(i)}{\rm{ln}}\chi_{k}-a_{k}\eta_{k}^{(i)}}\right],$
(34)

where the equality holds only when $\chi_{k}=\chi_{k}^{(i)}$.

Next, we optimize the power allocation to maximize the LB of the OF instead of
maximizing (26a) directly. In specific, the LB maximization problem to be
solved in the $i+1$-th iteration is formulated as

$\displaystyle\vspace{-0.8cm}\mathop{\max}\limits_{\\{\chi_{k}\\},\\{p_{k}^{p}\\},\\{p_{k}^{d}\\}}\;\;\;$
$\displaystyle\sum\nolimits_{k\in\mathcal{K}}{{\hat{w}}_{k}}^{(i)}\ln{\chi_{k}}$
(35a) $\displaystyle{\rm{s.t.}}\;\;\;$
$\displaystyle(\ref{djrdsfrefjgot}),(\ref{fhhthtyh}),(\ref{freocdsdi}),$ (35b)

where
${{\hat{w}}_{k}}^{(i)}={{\tilde{w}}_{k}}\hat{\rho}_{k}^{(i)}-a_{k}{{\tilde{w}}_{k}}\rho_{k}^{(i)}$
and the constant term in the OF is omitted.

Then, we can transform the above optimization into a GP problem as follows:

$\displaystyle\mathop{\max}\limits_{\\{\chi_{k}\\},\\{p_{k}^{p}\\},\\{p_{k}^{d}\\}}\;\;\;$
$\displaystyle\mathop{\prod}\nolimits_{k\in\mathcal{K}}\chi_{k}^{\hat{w}_{k}^{(i)}}$
(36a) $\displaystyle{\rm{s.t.}}\;\;\;$
$\displaystyle\sum\limits_{i\in\mathcal{K}\setminus
k}\\!\\!{{\alpha_{i}}{\alpha_{k}}{K}{\chi_{k}}p_{k}^{p}p_{i}^{d}}\\!\\!+\\!\\!\sum\limits_{i\in\mathcal{K}}{{\alpha_{i}}{\chi_{k}}p_{i}^{d}}\\!+\\!{\chi_{k}}{\alpha_{k}}{K}p_{k}^{p}\\!+\\!{\chi_{k}}\\!\leq\\!(M-1){K}\alpha_{k}^{2}p_{k}^{p}p_{k}^{d},\forall
k,$ (36b)
$\displaystyle{\chi_{k}}\geq{1\mathord{\left/{\vphantom{1{{f_{k}^{-1}}\left({\frac{{R_{k}^{{\rm{req}}}\ln
2}}{{1-\beta}}}\right)}}}\right.\kern-1.2pt}{{f_{k}^{-1}}\left({\frac{{R_{k}^{{\rm{req}}}\ln
2}}{{1-\beta}}}\right)}},\forall k,(\ref{freocdsdi}).$ (36c)

Although GP is not a convex optimization problem, it can be equivalently
transformed into a convex optimization problem by applying a logarithmic
change of variables. Therefore, the globally optimal solution of Problem (36)
can be obtained by using the interior-point method [27]. Some software
packages that can solve GP problem are MOSEK package and CVX [28].

Based on the above discussion, the iterative algorithm to solve Problem (26)
is given in Algorithm 1.

Algorithm 1 Iterative algorithm for solving Problem (26) for MRC

1: Initialize iteration number $i=1$, error tolerance $\xi$. Initialize a
feasible power allocation $\\{{p_{k}^{p}}^{(0)},{p_{k}^{d}}^{(0)},\forall
k\\}$, calculate
$\\{\chi_{k}^{(0)},\rho_{1}^{(0)},\rho_{2}^{(0)},{{\hat{w}}_{k}}^{(0)},\forall
k\\}$, and calculate the OF of Problem (26), denoted as ${\rm{Obj}}^{(0)}$.

2: With given
$\\{\chi_{k}^{(i-1)},\rho_{1}^{(i-1)},\rho_{2}^{(i-1)},{{\hat{w}}_{k}}^{(i-1)},\forall
k\\}$, solve Problem (36) by using the CVX package to obtain
$\\{p_{k}^{p(i)},p_{k}^{d(i)},\chi_{k}^{(i)},\forall k\\}$.

3: Update $\\{\rho_{k}^{(i)},\rho_{k}^{(i)},{{\hat{w}}_{k}}^{(i)},\forall
k\\}$;

4: Calculate the new OF ${\rm{Obj}}^{(i)}$. If
${{\left|{{\rm{Ob}}{{\rm{j}}^{(i)}}-{\rm{Ob}}{{\rm{j}}^{(i-1)}}}\right|}\mathord{\left/{\vphantom{{\left|{{\rm{Ob}}{{\rm{j}}^{(i)}}-{\rm{Ob}}{{\rm{j}}^{(i-1)}}}\right|}{{\rm{Ob}}{{\rm{j}}^{(i)}}}}}\right.\kern-1.2pt}{{\rm{Ob}}{{\rm{j}}^{(i)}}}}<\xi$,
terminate. Otherwise, set $i\leftarrow i+1$, go to step 2.

###  III-B Algorithm Analysis

####  III-B1 Initialization of Algorithm 1

As shown in Step 1 of Algorithm 1, one has to find a feasible initial power
allocation in order to make the algorithm work. Note that randomly selecting a
set of power allocation solutions that satisfy the per-device energy
constraints may not satisfy their minimum SINR requirements. Hence, one has to
carefully choose the initial power allocation. In the following, we provide
one alternative method to find the initial power allocation solutions.

Inspired by the user selection problem formulation in [29, 30], we construct
the following alternative optimization problem by introducing an auxiliary
variable $\varphi$:

$\displaystyle\mathop{\max}\limits_{\varphi,\\{p_{k}^{p}\\},\\{p_{k}^{d}\\}}\;\;\;$
$\displaystyle\varphi$ (37a) $\displaystyle{\rm{s.t.}}\;\;\;$
$\displaystyle{{\hat{\gamma}}_{k}}\geq{{\varphi\mathord{\left/{\vphantom{1{{f_{k}^{-1}}\left({\frac{{R_{k}^{{\rm{req}}}\ln
2}}{{1-\beta}}}\right)}}}\right.\kern-1.2pt}{{f_{k}^{-1}}\left({\frac{{R_{k}^{{\rm{req}}}\ln
2}}{{1-\beta}}}\right)}}},\forall k,(\ref{freocdsdi}).$ (37b)

Obviously, Problem (37) is always feasible since at least
$\\{\varphi=0,p_{k}^{p}=0,p_{k}^{d}=0,\forall k\\}$ is a feasible solution. It
can be readily verified that the original Problem (24) is feasible if the
optimal $\varphi\geq 1$, and the output power allocation can be adopted as the
initial input for Algorithm 1. In this paper, we assume that Problem (24) is
always feasible and the optimal $\varphi$ in Problem (37) is always no smaller
than one. Problem (37) can also be transformed into a GP problem, where the
globally optimal solution can be obtained. The details are omitted here due to
the limited space.

####  III-B2 Convergence Analysis

In this part, we analyze the convergence of Algorithm 1. In the following, we
show that ${\rm{Obj}}^{(i)}\leq{\rm{Obj}}^{(i+1)}$.

Since $\\{\chi_{k}^{(i+1)},\forall k\\}$ is the optimal solution of Problem
(35) in the $i+1$-th iteration, we have

$\begin{array}[]{l}\quad\sum\nolimits_{k\in\mathcal{K}}{{{\tilde{w}}_{k}}}\left[{\hat{\rho}_{k}^{(i)}\ln\left({\chi_{k}^{(i+1)}}\right)+\hat{\eta}_{k}^{(i)}-{a_{k}}\rho_{k}^{(i)}{\rm{ln}}(\chi_{k}^{(i+1)})-{a_{k}}\eta_{k}^{(i)}}\right]\\\
\geq\sum\nolimits_{k\in\mathcal{K}}{{{\tilde{w}}_{k}}}\left[{\hat{\rho}_{k}^{(i)}\ln\left({\chi_{k}^{(i)}}\right)+\hat{\eta}_{k}^{(i)}-{a_{k}}\rho_{k}^{(i)}{\rm{ln}}(\chi_{k}^{(i)})-{a_{k}}\eta_{k}^{(i)}}\right]={\rm{Obj}}^{(i)}.\end{array}$
(38)

By using inequality (34) with $\chi_{k}=\chi_{k}^{(i+1)}$, we have

$\begin{array}[]{l}\quad\sum\nolimits_{k\in\mathcal{K}}{{\tilde{w}}_{k}}\left[{\ln(1+\chi_{k}^{(i+1)})-{a_{k}}G(\chi_{k}^{(i+1)})}\right]\\\
\geq\sum\nolimits_{k\in\mathcal{K}}{{\tilde{w}}_{k}}\left[{\hat{\rho}_{k}^{(i)}\ln(\chi_{k}^{(i+1)})+\hat{\eta}_{k}^{(i)}-{a_{k}}\rho_{k}^{(i)}{\rm{ln}}(\chi_{k}^{(i+1)})-{a_{k}}\eta_{k}^{(i)}}\right].\end{array}$
(39)

By combining (38) and (39), we have

${\rm{Ob}}{{\rm{j}}^{(i+1)}}\geq\sum\nolimits_{k\in\mathcal{K}}{{{\tilde{w}}_{k}}}\left[{\hat{\rho}_{k}^{(i)}\ln(\chi_{k}^{(i+1)})+\hat{\eta}_{k}^{(i)}-{a_{k}}\rho_{k}^{(i)}{\rm{ln}}(\chi_{k}^{(i+1)})-{a_{k}}\eta_{k}^{(i)}}\right]\geq{\rm{Ob}}{{\rm{j}}^{(i)}}.$
(40)

In addition, since each device has its own energy constraint, the OF value of
Problem (24) has upper bound. As a result, Algorithm 1 is guaranteed to
converge.

####  III-B3 Solution Analysis

Since the original problem (24) is non-convex, the globally optimal solution
cannot be obtained in general. However, by using the similar proof as in
Appendix B in [30], we can prove that Algorithm 1 can converge to the Karush-
Kuhn-Tucker (KKT) point of Problem (24). The converged solution only depends
on the initial input of Algorithm 1. Simulation results show that the
algorithm almost converges to the same solution with different initial input
solutions.

####  III-B4 Complexity Analysis

The main complexity mainly lies in solving a GP problem in each iteration. In
[31], the authors claimed that the GP problem can be efficiently solved by
using the standard interior point methods with a worst-case polynomial-time
complexity. The upper bound of the total number of Newton steps in the
interior point method does not depend on the number of variables, or the
number of constraints. The derived upper bound shows that the barrier method
converges linearly. By carefully choosing the parameters, the bound on the
number of Newton steps can grow as $\sqrt{m}$ instead of $m$, where $m$ is the
number of constraints [32]. In addition, simulation results show that
Algorithm 1 converges rapidly, which means that Algorithm 1 can converge to a
local optimal solution with a polynomial time complexity.

##  IV Weighted Sum Data Rate for ZF

In this section, we aim to deal with weighted sum rate maximization problem
for ZF. Due to different expressions of the SINR of MRC and ZF, some
derivations in MRC cannot be directly applied in the ZF case. In the
following, we develop an efficient algorithm to solve Problem (26) for the ZF
case.

###  IV-A Algorithm Design

By substituting (5) and (6) into (23), the SINR in the ZF case can be
reformulated as

${{\hat{\gamma}}_{k}}=\frac{{(M-K)\alpha_{k}^{2}Kp_{k}^{p}p_{k}^{d}\prod\nolimits_{i\in\mathcal{K}}{(1+{\alpha_{i}}Kp_{i}^{p})}}}{{\left({1+{\alpha_{k}}Kp_{k}^{p}}\right)\left({\sum\nolimits_{i\in\mathcal{K}}{p_{i}^{d}{\alpha_{i}}\prod\nolimits_{j\neq
i}{(1+{\alpha_{j}}Kp_{j}^{p})}+\prod\nolimits_{i\in\mathcal{K}}{(1+{\alpha_{i}}Kp_{i}^{p})}}}\right)}},\forall
k.$ (41)

Unfortunately, since the nominator of ${{\hat{\gamma}}_{k}}$ in (41) is a
posynomial function, the SINR constraint in (26b) cannot be transformed into
the format as that in (36b), in which the right hand side is a monomial
function. Hence, the problem cannot be directly transformed into a GP problem.
To solve this issue, we introduce Theorem 3 as follows.

_Theorem 3_ : For any given vector ${\breve{\bf
x}}=\\{\breve{x}_{1},\cdots,\breve{x}_{K}\\}$ with $\breve{x}_{i}\geq
0,\forall n$, function
$W({\bf{x}})=\prod\nolimits_{i\in\mathcal{K}}{(1+{x_{i}})}$ is lower bounded
by

$\vspace{-0.1cm}W({\bf{x}})=\prod\nolimits_{i\in\mathcal{K}}{(1+{x_{i}})}\geq\lambda\prod\nolimits_{i\in\mathcal{K}}{x_{i}^{{\tau_{i}}}}\buildrel\Delta\over{=}Y({\bf{x}})$
(42)

where $\lambda$ and $\tau_{i},\forall i$ are given by

$\lambda=\frac{{\prod\nolimits_{i\in\mathcal{K}}{(1+{{\breve{x}}_{i}})}}}{\prod\nolimits_{i\in\mathcal{K}}{\breve{x}_{i}^{{\tau_{i}}}}},{\tau_{i}}=\frac{{{{\breve{x}}_{i}}}}{{1+{{\breve{x}}_{i}}}},\forall
i,$ (43)

and ${\bf{x}}$ is given by ${\bf{x}}=\\{x_{1},\cdots,x_{K}\\}$.

In addition, we have:

$\vspace{-0.45cm}W({\breve{\bf x}})=Y({\breve{\bf x}}),\nabla W({\breve{\bf
x}})=\nabla Y({\breve{\bf x}}).$ (44)

where $\nabla W(\bf{x})$ and $\nabla Y(\bf{x})$ denote the gradient of
function $W(\cdot)$ and $Y(\cdot)$ w.r.t. $\bf{x}$, respectively.

_Proof_ : Please refer to Appendix E.

Based on Theorem 3, we replace the polynomial functions in the nominator of
${{\hat{\gamma}}_{k}}$ in (41) with their best local monomial approximations
by employing Theorem 3 222The best local monomial approximations means that
the approximation function should satisfy three conditions as specified in
Section IV-A of [33].. In specific, we denote
${\bf{P}}^{(n)}=\\{p_{k}^{p(n)},p_{k}^{d(n)},\forall k\\}$ as the power
allocation in the $n$-th iteration, and the corresponding $\chi_{k}$ is given
by $\chi_{k}^{(n)}$. Then, in the $n+1$-th iteration, we approximate the term
${\prod\nolimits_{i\in\mathcal{K}}{(1+{\alpha_{i}}Kp_{i}^{p})}}$ in the
nominator of (41) by its best local monomial approximations, which is given by
function $Y({\bf{x}})$ in Theorem 3 with $x_{i}={\alpha_{i}}Kp_{i}^{p}$:

$\prod\nolimits_{i\in\mathcal{K}}{(1+{\alpha_{i}}Kp_{i}^{p})}\geq{\lambda^{(n)}}\prod\nolimits_{i\in\mathcal{K}}{{{\left({{\alpha_{i}}Kp_{i}^{p}}\right)}^{{\tau_{i}^{(n)}}}}},$
(45)

where ${\lambda^{(n)}}$ and ${\tau_{i}^{(n)}},\forall i$ are given in (43)
with ${\tilde{x}}_{i}={\alpha_{i}}Kp_{i}^{p(n)},\forall i$. Then, we focus on
the following constraint instead of the original SINR constraint in (26b):

$\displaystyle{\chi_{k}}\left({1+{\alpha_{k}}Kp_{k}^{p}}\right)\left({\sum\nolimits_{i\in\mathcal{K}}{p_{i}^{d}{\alpha_{i}}\prod\nolimits_{j\neq
i}{(1+{\alpha_{j}}Kp_{j}^{p})}+\prod\nolimits_{i\in\mathcal{K}}{(1+{\alpha_{i}}Kp_{i}^{p})}}}\right)$
$\displaystyle\leq$
$\displaystyle(M-K)\alpha_{k}^{2}Kp_{k}^{p}p_{k}^{d}{\lambda^{(n)}}\prod\nolimits_{i\in\mathcal{K}}{{{\left({{\alpha_{i}}Kp_{i}^{p}}\right)}^{{\tau_{i}}^{(n)}}}}.$
(46)

Note that the left hand side (LHS) of (IV-A) is a posynomial function, and the
right hand side (RHS) becomes a monomial function. In addition, the LHS is no
larger than RHS. Hence, constraint (IV-A) satisfies the conditions for a
problem to be a GP problem [27].

Then, by using the same method as in the MRC case to deal with the OF, in the
$n+1$-th iteration, we aim to solve the following GP problem:

$\displaystyle\vspace{-0.2cm}\mathop{\max}\limits_{\\{\chi_{k}\\},\\{p_{k}^{p}\\},\\{p_{k}^{d}\\}}\;\;\;$
$\displaystyle\mathop{\prod}\nolimits_{k\in\mathcal{K}}\chi_{k}^{\hat{w}_{k}^{(n)}}$
(47a) $\displaystyle{\rm{s.t.}}\;\;\;$
$\displaystyle(\ref{freocdsdi}),(\ref{ilhiuytrew})$ (47b)
$\displaystyle{{\hat{\gamma}}_{k}}\geq{1\mathord{\left/{\vphantom{1{{f_{k}^{-1}}\left({\frac{{R_{k}^{{\rm{req}}}\ln
2}}{{1-\beta}}}\right)}}}\right.\kern-1.2pt}{{f_{k}^{-1}}\left({\frac{{R_{k}^{{\rm{req}}}\ln
2}}{{1-\beta}}}\right)}},\forall k,$ (47c)

where the parameters ${\hat{w}_{k}^{(n)}}$’s are the same as those in the MRC
case. This problem can be efficiently solved by using CVX [28].

Based on the above discussion, the iterative algorithm to solve Problem (26)
for the case of ZF is provided in Algorithm 2.

Algorithm 2 Iterative algorithm for solving Problem (26) for ZF

1: Initialize iteration number $n=1$, error tolerance $\xi$. Initialize a
feasible power allocation $\\{p_{k}^{p(0)},p_{k}^{d(0)},\forall k\\}$,
calculate
$\\{\chi_{k}^{(0)},\rho_{1}^{(0)},\rho_{2}^{(0)},{{\hat{w}}_{k}}^{(0)},\lambda^{(0)},\tau_{k}^{(0)},\forall
k\\}$, and calculate the OF of Problem (26), denoted as ${\rm{Obj}}^{(0)}$.

2: With given
$\\{\chi_{k}^{(n-1)},\rho_{1}^{(n-1)},\rho_{2}^{(n-1)},{{\hat{w}}_{k}}^{(n-1)},\lambda^{(n-1)},\tau_{k}^{(n-1)}\forall
k\\}$, solve Problem (47) by using the CVX package to obtain
$\\{p_{k}^{p(n)},p_{k}^{d(n)},\chi_{k}^{(n)},\forall k\\}$.

3: Update
$\\{\rho_{1}^{(n)},\rho_{2}^{(n)},{{\hat{w}}_{k}}^{(n)},\lambda^{(n)},\tau_{k}^{(n)},\forall
k\\}$;

4: Calculate new OF ${\rm{Obj}}^{(n)}$. If
${{\left|{{\rm{Ob}}{{\rm{j}}^{(n)}}-{\rm{Ob}}{{\rm{j}}^{(n-1)}}}\right|}\mathord{\left/{\vphantom{{\left|{{\rm{Ob}}{{\rm{j}}^{(n)}}-{\rm{Ob}}{{\rm{j}}^{(n-1)}}}\right|}{{\rm{Ob}}{{\rm{j}}^{(n)}}}}}\right.\kern-1.2pt}{{\rm{Ob}}{{\rm{j}}^{(n)}}}}<\xi$,
terminate. Otherwise, set $n\leftarrow n+1$, go to step 2.

###  IV-B Algorithm Analysis

Algorithm 2 can be analyzed similar to Algorithm 1 except the convergence
analysis since the problems to be solved in each iteration of Algorithm 2 do
not have the same set of constraints.

In the following, we prove that the solution obtained in the $n$-th iteration
is also feasible for the problem to be solved in the $n+1$-th iteration. We
only need to check constraint (IV-A) since the other two constraints are the
same in each iteration.

Let us denote $\\{\chi_{k}^{(n)},p_{k}^{p(n)},p_{k}^{d(n)},\forall k\\}$ as
the optimal solution in the $n$-th iteration. Then, it is also a feasible
solution, and we have

$\displaystyle\chi_{k}^{(n)}\left({1+{\alpha_{k}}Kp_{k}^{p(n)}}\right)\left({\sum\nolimits_{i\in\mathcal{K}}{p_{i}^{d(n)}{\alpha_{i}}\prod\nolimits_{j\neq
i}{(1+{\alpha_{j}}Kp_{j}^{p(n)})}+\prod\nolimits_{i\in\mathcal{K}}{(1+{\alpha_{i}}Kp_{i}^{p(n)})}}}\right)$
$\displaystyle\leq$
$\displaystyle(M-K)\alpha_{k}^{2}Kp_{k}^{p(n)}p_{k}^{d(n)}{\lambda^{(n-1)}}\prod\nolimits_{i\in\mathcal{K}}{{{\left({{\alpha_{i}}Kp_{i}^{p(n)}}\right)}^{\tau_{i}^{(n-1)}}}}.$
(48)

By using (45), we have

$\prod\nolimits_{i\in\mathcal{K}}{(1+{\alpha_{i}}Kp_{i}^{p(n)})}\geq{\lambda^{(n-1)}}\prod\nolimits_{i\in\mathcal{K}}{{{\left({{\alpha_{i}}Kp_{i}^{p(n)}}\right)}^{\tau_{i}^{(n-1)}}}}.$
(49)

In addition, by using (44) in Theorem 3 and (49), we have

${\lambda^{(n)}}\prod\nolimits_{i\in\mathcal{K}}{{{\left({{\alpha_{i}}Kp_{i}^{p(n)}}\right)}^{\tau_{i}^{(n)}}}}=\prod\nolimits_{i\in\mathcal{K}}{(1+{\alpha_{i}}Kp_{i}^{p(n)})}\geq{\lambda^{(n-1)}}\prod\nolimits_{i\in\mathcal{K}}{{{\left({{\alpha_{i}}Kp_{i}^{p(n)}}\right)}^{\tau_{i}^{(n-1)}}}}.$
(50)

Finally, by combining (IV-A) and (50), we have

$\displaystyle\chi_{k}^{(n)}\left({1+{\alpha_{k}}Kp_{k}^{p(n)}}\right)\left({\sum\nolimits_{i\in\mathcal{K}}{p_{i}^{d(n)}{\alpha_{i}}\prod\nolimits_{j\neq
i}{(1+{\alpha_{j}}Kp_{j}^{p(n)})}+\prod\nolimits_{i\in\mathcal{K}}{(1+{\alpha_{i}}Kp_{i}^{p(n)})}}}\right)$
$\displaystyle\leq$
$\displaystyle(M-K)\alpha_{k}^{2}Kp_{k}^{p(n)}p_{k}^{d(n)}{\lambda^{(n)}}\prod\nolimits_{i\in\mathcal{K}}{{{\left({{\alpha_{i}}Kp_{i}^{p(n)}}\right)}^{\tau_{i}^{(n)}}}}.$
(51)

Hence, $\\{\chi_{k}^{(n)},p_{k}^{p(n)},p_{k}^{d(n)},\forall k\\}$ is also a
feasible solution in the $n+1$-th iteration. Then, by using the similar proof
as in the case of MRC, we can also prove that Algorithm 2 is guaranteed to
converge. By using a similar method, we also prove that this algorithm will
converge to a feasible solution of Problem (26).

##  V Simulation Results

In this section, we provide simulation results to demonstrate the
effectiveness of our proposed algorithms for industrial automation systems.
The channel path loss is modeled as $PL=35.3+37.6{\rm{log}}_{10}d$ (dB) [34],
and the small-scale fading is modeled as Rayleigh fading with zero mean and
unit variance. Unless otherwise specified, the simulation parameters are set
as follows: number of transmit antennas of $M=100$, number of devices of
$K=10$, channel bandwidth of $B=0.2$ MHz, noise power spectral density of -174
dBm/Hz, decoding error probability of ${\varepsilon_{k}}=10^{-9},\forall k$,
number of transmit antennas of $M=100$, and maximum transmission duration of
$0.5$ ms. The maximum blocklength is calculated as $L=BT=100$. The other
parameters are specified in each figure. The energy constraint for each device
is assumed to be equal, i.e., $E_{k}=E,\forall k$, and each device has the
same data rate targets, $R_{k}^{\rm{req}}=R^{\rm{req}},\forall k$. The weights
for each device are uniformly generated within $[0,1]$.

###  V-A Tightness of the Date Rate LB

In Fig. 6 and Fig. 6, we investigates the tightness of the LB derived for the
cases of MRC and ZF, respectively. The simulation results are obtained through
the Monte-Carlo simulation by averaging over 5000 random channel generations.
It is observed that the data rate LB is tight for both cases of MRC and ZF for
any number of transmit antennas. Interestingly, the curves for the ZF case are
almost overlapped with each other. This verifies that the data rate LBs
derived in Theorem 1 and Theorem 2 are suitable for optimization instead of
directly optimizing the complicated expectation expression.

Figure 5: Tightness of the derived data rate LB for the MRC case.

Figure 6: Tightness of the derived data rate LB for the ZF case.

###  V-B Convergence Behaviour of the Proposed Algorithms

Fig. 8 and Fig. 8 investigate the impact of the energy limit at each device on
the performance of the proposed Algorithm 1 for the MRC case and Algorithm 2
for the ZF case, respectively. From these two figures, we can observe that
both algorithms converge rapidly for various energy limits and only 2 or 3
iterations are sufficient for the algorithms to converge. This demonstrates
the low complexity of the proposed algorithms.

Figure 7: Convergence behaviour of the proposed algorithm for the MRC case.

Figure 8: Convergence behaviour of the proposed algorithm for the ZF case.

###  V-C Performance Comparison

In this subsection, we compare the proposed algorithms with the following
algorithms:

  * •

Upper bound: In this method, Shannon capacity is adopted for optimization in
Problem (26). In other words, the penalty term $a_{k}G({\chi_{k}})$ is set to
zero in both the OF and rate constraint (24b), i.e.,
$a_{k}G({\chi_{k}})=0,\forall k$. This method provides the upper bound of the
average weighted sum rate performance in the IIoT networks.

  * •

Conventional alg.: As in [35], the solution obtained from the upper bound is
applied in (20) and (22) for the cases of MRC and ZF respectively by
considering the penalty term $a_{k}G({\chi_{k}})$. That means that the upper
bound is used for obtaining solutions, but the achievable data rate under
finite blocklength is used for performance evaluation.

  * •

Fixed pilot power alg.: In this scheme, we only optimize the payload power
$p_{k}^{d}$ while fixing the pilot power as $p_{k}^{p}=E/L$. This algorithm is
provided to show the benefits of jointly optimizing the pilot power and
payload power.

The following results are obtained by averaging over 100 Monte-Carlo
simulations where in each snapshot the devices are randomly generated in the
cell. For each snapshot, if the device’s achievable data rate cannot achieve
its rate targets, we set the corresponding data rate to zero. The proposed
algorithm is denoted as ‘Proposed alg.’ in the following figures.

In Fig. 10 and Fig. 10, we show the average weighted sum rate versus the
energy limit at each device, $E$. The rate targets for MRC and ZF are set as
$R^{\rm{req}}=1$ bit/s/Hz and $R^{\rm{req}}=4$ bit/s/Hz, respectively. As seen
from these figures, the system performance increases with the available energy
at each device since the SINR at each device is increased. Note that the
weighted sum rate achieved by some algorithms may approach zero, which means
the power allocation solution is infeasible. As expected, the upper bound has
the best performance since the penalty terms are not considered. Furthermore,
the ‘conventional alg.’ based on Shannon capacity formula has higher
probability to violate the data rate requirement, especially for samll $E$.
Therefore, Shannon capacity cannot be employed for the transmission design of
URLLC for industrial applications, in particular when the energy limit is
small as the QoS target cannot be guaranteed. However, for the large value of
$E$, the SINR value for each device is very high and then the penalty term
$a_{k}G({\chi_{k}})$ can be ignored when comparing with the first term of
${\ln}(1+{\chi_{k}})$. Hence, the ‘conventional alg.’ will achieve similar
performance as that of the proposed algorithm. As shown in Fig. 10 and Fig.
10, by jointly optimizing the payload power and pilot power, the proposed
algorithm is superior over the ‘Fixed pilot power alg.’, which only optimizes
the payload power. The performance gain is obvious when the energy limit is
low, especially for the case of ZF. This means the need of jointly optimizing
pilot and payload power at low energy limit. This can be explained as follows.
When the energy limit is low, the system performance is limited by channel
estimation procedure. Through the joint power allocation, some power can be
borrowed from that for data transmission to enhance the channel estimation
accuracy, and thus increases the weighted sum rate. However, in the high
energy $E$, the accuracy of channel estimation is already enough, and
additional joint power control bring marginal performance gain. Another
interesting observation is that the performance gain of the proposed algorithm
over ‘Fixed pilot power alg.’ is much smaller for the MRC case than the ZF
case. This may be due to the fact that the CSI is more important when using it
for removing multiple device interference in the ZF case.

Figure 9: Average weighted sum rate vs. maximum energy limit
$10{\rm{log}}_{10}E$ for various schemes for the case of MRC.

Figure 10: Average weighted sum rate vs. maximum energy limit
$10{\rm{log}}_{10}E$ for various schemes for the case of ZF.

Figure 11: Average weighted sum rate vs. the number of devices for various
schemes for the case of MRC.

Figure 12: Average weighted sum rate vs. the number of devices for various
schemes for the case of ZF.

Fig. 12 and Fig. 12 show the average weighted sum rate versus the number of
devices for various schemes. The rate targets for MRC and ZF are set as
$R^{\rm{req}}=1$ bit/s/Hz and $R^{\rm{req}}=2$ bit/s/Hz, respectively. The
energy limit for MRC and ZF are set as $E=2$ and $E=1$, respectively. It is
observed from Fig. 12 and Fig. 12 that the average weighted sum rate achieved
by all the algorithms except the ‘Conventional alg.’ increases with the number
of devices since these schemes fully exploit the multi-device diversity. For
the MRC case, the performance of the ‘Conventional alg.’ decreases with the
number of devices. The main reason is that the design is based on Shannon
capacity formula, which does not take into consideration the effect of short
blocklength on the achievable data rate. As a result, the probability that the
achieved data rate violates the rate requirement will increase with the number
of devices. It is interesting to observe that ‘Fixed pilot power alg.’ has
similar performance as the proposed algorithm for the case of MF. However, for
the case of ZF, the proposed algorithm significantly outperforms ‘Fixed pilot
power alg.’, and the performance gain is increasing with the number of
devices. This again reveals the importance of joint power optimization in the
case of ZF. For the ZF case, the weighted sum rate achieved by ‘Conventional
alg.’ first increases with the number of devices and then decreases with it.

Finally, we study the effect of blocklength on the weighted sum rate
performance in Fig. 14 and Fig. 14 for the cases of MRC and ZF, respectively.
The rate targets for MRC and ZF are set as $R^{\rm{req}}=2$ bit/s/Hz and
$R^{\rm{req}}=4$ bit/s/Hz, respectively. As expected, the system performance
increases with the channel blocklength since more time/frequency resource can
be exploited for transmission. Some interesting observations can be found in
these figures. First, when the blocklength is small, there is a significant
performance gap between the proposed algorithm and the upper bound. However,
this gap starts to reduce with the increase of the blocklength. This
demonstrates that the blocklength has significant impact for IIoT devices with
URLLC requirements. In addition, the proposed algorithm outperforms the ‘Fixed
pilot power alg.’ algorithm for both cases of MRC and ZF. This may be due to
the fact that higher rate target is imposed in this example than that in Fig.
12.

Figure 13: Average weighted sum rate vs. blocklength for various schemes for
the case of MRC.

Figure 14: Average weighted sum rate vs. blocklength for various schemes for
the case of ZF.

##  VI Conclusions

In this paper, we studied the resource allocation for uplink massive MIMO
systems to support critical IIoT operating under finite channel blocklength,
where multiple robots and/or actuators transmit URLLC signals to the central
controller simultaneously. We first derived the closed-form data rate LB with
imperfect CSI for both MRC and ZF receivers under the short packet
transmission. We then formulate the weighted sum rate maximization to jointly
optimize the pilot and payload power allocation while considering their
energy, minimum data rate and decoding error probability requirements. This
optimization problem is non-convex, and we proposed novel low-complexity
iterative algorithms to solve it. Simulation results demonstrate that the
algorithm converges rapidly, and outperforms the existing benchmark
algorithms, especially the algorithm based on conventional Shannon capacity.
This reveals the importance of adopting the achievable data expression for
finite channel blocklength.

##  Appendix A Proof of Theorem 1

We first consider the MRC beamforming. The proof follows the similar steps as
those in Appendix A of [36] for perfect CSI. Denote $\gamma_{k}$ as the
instantaneous SINR value when using MRC. By substituting
${{\bf{a}}_{k}}={{{\hat{\bf h}}}_{k}}$ into (12), we have

$\gamma_{k}=\frac{{p_{k}^{d}{{\left\|{{{{\hat{\bf
h}}}_{k}}}\right\|}^{4}}}}{{\sum\nolimits_{i\in\mathcal{K}\setminus
k}{p_{i}^{d}{{\left|{{\hat{\bf h}}_{k}^{H}{{{\hat{\bf
h}}}_{i}}}\right|}^{2}}}+\sum\nolimits_{i\in\mathcal{K}}{p_{i}^{d}{{\left|{{\hat{\bf
h}}_{k}^{H}{{{\tilde{\bf h}}}_{i}}}\right|}^{2}}}+{{\left\|{{{{\hat{\bf
h}}}_{k}}}\right\|}^{2}}}}.$ (52)

Then, ${{\mathbb{E}}\left\\{{\frac{1}{{{\gamma_{k}}}}}\right\\}}$ can be
expressed as

${\mathbb{E}}\left\\{{\frac{1}{{\gamma_{k}}}}\right\\}={\mathbb{E}}\left\\{{\frac{{\sum\nolimits_{i\in\mathcal{K}\setminus
k}{p_{i}^{d}{{\left|{{u_{k,i}}}\right|}^{2}}}+\sum\nolimits_{i\in\mathcal{K}}{p_{i}^{d}{{\left|{{v_{k,i}}}\right|}^{2}}}+1}}{{p_{k}^{d}{{\left\|{{{{\bf{\hat{h}}}}_{k}}}\right\|}^{2}}}}}\right\\}$
(53)

where ${u_{k,i}}=\frac{{{\hat{\bf
h}}_{k}^{H}{{{\bf{\hat{h}}}}_{i}}}}{{{{\left\|{{{{\bf{\hat{h}}}}_{k}}}\right\|}}}}$
and ${v_{k,i}}=\frac{{{\hat{\bf h}}_{k}^{H}{{{\tilde{\bf
h}}}_{i}}}}{{{{\left\|{{{{\hat{\bf h}}}_{k}}}\right\|}}}}$. Conditioned on
${\hat{\bf h}}_{k}$, ${u_{k,i}}$ and ${v_{k,i}}$ are Gaussian random variables
with zero mean and variance equal to $\sigma_{i}$ and $\delta_{i}$,
respectively. In addition, ${u_{k,i}}$ and ${v_{k,i}}$ are independent of
${{\hat{\bf h}}}_{k}$. Then, we have

${\mathbb{E}}\left\\{{\frac{1}{{\gamma_{k}}}}\right\\}=\left({\sum\nolimits_{i\in\mathcal{K}\setminus
k}{p_{i}^{d}{\sigma_{i}}}+\sum\nolimits_{i\in\mathcal{K}}{p_{i}^{d}{\delta_{i}}}+1}\right){\mathbb{E}}\left\\{{\frac{1}{{p_{k}^{d}{{\left\|{{{{\hat{\mathbf{h}}}}_{k}}}\right\|}^{2}}}}}\right\\}.$
(54)

By using the identity [37]

$\vspace{-0.2cm}{\mathbb{E}}\left\\{{{\rm{tr}}\left({{{\bf{W}}^{-1}}}\right)}\right\\}=\frac{m}{{n-m}}$
(55)

where ${\bf{W}}\sim{{\cal W}_{m}}\left({n,{{\bf{I}}_{n}}}\right)$ is an
$m\times m$ central complex Wishart matrix with $n$ ($n-m$) degrees of
freedom. Then, we have

${\mathbb{E}}\left\\{{\frac{1}{{p_{k}^{d}{{\left\|{{{{\hat{\mathbf{h}}}}_{k}}}\right\|}^{2}}}}}\right\\}=\frac{1}{{p_{k}^{d}(M-1){\sigma_{k}}}},\
{\rm{for}}\ M\geq 2.$ (56)

By substituting (56) into (54), we can obtain ${{\hat{\gamma}}_{k}}$ in (21).

##  Appendix B Proof of Theorem 2

By using ZF, we have ${{\bf{A}}^{H}}{\hat{\bf H}}={{\bf{I}}_{K}}$. Thus,
${{\bf{a}}_{k}}{{{\hat{\bf h}}}_{i}}$ is equal to one if $k=i$; otherwise, it
is equal to zero. Then, the instantaneous SINR for ZF can be rewritten as

$\gamma_{k}=\frac{{p_{k}^{d}}}{{\sum\nolimits_{i\in\mathcal{K}}{p_{i}^{d}{{\left|{{\bf{a}}_{k}^{H}{{{\tilde{\bf
h}}}_{i}}}\right|}^{2}}}+{{\left\|{{{\bf{a}}_{k}}}\right\|}^{2}}}}.$ (57)

Define
${\bf{\Lambda}}={\rm{diag}}\left\\{{\sqrt{{\sigma_{1}}},\cdots,\sqrt{{\sigma_{K}}}}\right\\}$
and ${\bf{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\smile$}}\over{H}}}}=\left[{\frac{{{{{\bf{\hat{h}}}}_{1}}}}{{\sqrt{{\sigma_{1}}}}},\cdots,\frac{{{{{\bf{\hat{h}}}}_{K}}}}{{\sqrt{{\sigma_{K}}}}}}\right]$.
Then, we can express ${\bf{\hat{H}}}$ as
${\bf{\hat{H}}}={\bf{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\smile$}}\over{H}}\Lambda}}$. The columns of
${\bf{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\smile$}}\over{H}}}}$ are independent of each
other, and each column follows the distribution of
$\mathcal{CN}(\mathbf{0},\mathbf{I})$. Then,
${{\mathbb{E}}\left\\{{\frac{1}{{{\gamma_{k}}}}}\right\\}}$ can be expressed
as

$\displaystyle{\mathbb{E}}\left\\{{\frac{1}{{\gamma_{k}}}}\right\\}$
$\displaystyle=$
$\displaystyle\left.{{\left({\sum\nolimits_{i\in\mathcal{K}}{p_{i}^{d}{\delta_{i}}+1}}\right)\mathbb{E}\left\\{{{{\left[{{{\left({{{{\bf{\hat{H}}}}^{H}}{\bf{\hat{H}}}}\right)}^{-1}}}\right]}_{k,k}}}\right\\}}}\middle/{{p_{k}^{d}}}\right.$
(58) $\displaystyle=$
$\displaystyle\left.{{\left({\sum\nolimits_{i\in\mathcal{K}}{p_{i}^{d}{\delta_{i}}+1}}\right){\mathbb{E}}\left\\{{{{\left[{{{\left({{\bf{\Lambda}}{{{\bf{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\smile$}}\over{H}}}}}^{H}}{\bf{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\smile$}}\over{H}}\Lambda}}}\right)}^{-1}}}\right]}_{k,k}}}\right\\}}}\middle/{{p_{k}^{d}}}\right.$
(59) $\displaystyle=$
$\displaystyle\left.{{\left({\sum\nolimits_{i\in\mathcal{K}}{p_{i}^{d}{\delta_{i}}+1}}\right){\mathbb{E}}\left\\{{{{\left[{{{\left({{{{\bf{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\smile$}}\over{H}}}}}^{H}}{\bf{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\smile$}}\over{H}}}}}\right)}^{-1}}}\right]}_{k,k}}}\right\\}}}\middle/{{{\sigma_{k}}p_{k}^{d}}}\right.$
(60) $\displaystyle=$
$\displaystyle\left.{{\left({\sum\nolimits_{i\in\mathcal{K}}{p_{i}^{d}{\delta_{i}}+1}}\right){\mathbb{E}}\left\\{{{\rm{tr}}\left[{{{\left({{{{\bf{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\smile$}}\over{H}}}}}^{H}}{\bf{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\smile$}}\over{H}}}}}\right)}^{-1}}}\right]}\right\\}}}\middle/{{K{\sigma_{k}}p_{k}^{d}}}\right.$
(61) $\displaystyle=$
$\displaystyle\left.{{\left({\sum\nolimits_{i\in\mathcal{K}}{p_{i}^{d}{\delta_{i}}+1}}\right)}}\middle/{{(M-K){\sigma_{k}}p_{k}^{d}}}\right.$
(62)

where (62) is obtained by using (55).

##  Appendix C Proof of Lemma 2

For notational simplicity, we define function
$T(x)\buildrel\Delta\over{=}1-\frac{1}{{{{(1+{x})}^{2}}}}$. The second
derivative of $T(x)$ w.r.t. $x$ is given by

$\vspace{-0.2cm}T^{\prime\prime}(x)=\frac{{-6}}{{{{\left({{x}+1}\right)}^{4}}}}\leq
0.$ (63)

Next, we prove that $G(x)$ is a concave function of $x$. Since $T(x)$ is a
concave function w.r.t. $x$, we have

$\vspace{-0.3cm}\theta T(\hat{x})+(1-\theta)T(\tilde{x})\leq
T(\theta\hat{x}+(1-\theta)\tilde{x}),$ (64)

for $0\leq\theta\leq 1$, where $\hat{x}$ and $\tilde{x}$ are two different
non-negative values. On the other hand, $\sqrt{x}$ is a concave function
w.r.t. $x$. Then, for $0\leq\theta\leq 1$ we have

$\vspace{-0.2cm}\theta\sqrt{T(\hat{x})}+(1-\theta)\sqrt{T(\tilde{x})}\leq\sqrt{\theta
T(\hat{x})+(1-\theta)T(\tilde{x})}.$ (65)

By combining (65) with (64), we have

$\theta\sqrt{T(\hat{x})}+(1-\theta)\sqrt{T(\tilde{x})}\leq\sqrt{T(\theta\hat{x}+(1-\theta)\tilde{x})},$
(66)

which is equivalent to $\theta G(\hat{x})+(1-\theta)G(\tilde{x})\leq
G(\theta\hat{x}+(1-\theta)\tilde{x})$. Hence, $G(x)$ is a concave function
w.r.t. $x$, which completes the proof.

##  Appendix D Proof of Lemma 3

The equalities in (30) can be readily proved by substituting the expressions
of $\rho$ and $\eta$ in (28) and (29) into (27). Next, we focus on the proof
of Inequality (27).

Define function $H(x)\buildrel\Delta\over{=}F(x)-G(x)$, and thus
$H(\tilde{x})=0$. Then, with some simple manipulations, the first-order
derivative of function $H(x)$ w.r.t. $x$ can be calculated as

$H^{\prime}(x)=\frac{{\tilde{x}\sqrt{{x^{2}}+2x}{{\left({x+1}\right)}^{2}}-x\sqrt{{{\tilde{x}}^{2}}+2\tilde{x}}{{\left({\tilde{x}+1}\right)}^{2}}}}{{\sqrt{{{\tilde{x}}^{2}}+2\tilde{x}}{{\left({\tilde{x}+1}\right)}^{2}}x\sqrt{{x^{2}}+2x}{{\left({x+1}\right)}^{2}}}}.$
(67)

Since both $x$ and $\tilde{x}$ are positive values, the sign of
$H^{\prime}(x)$ only depends on the nominator of $H^{\prime}(x)$. Then, denote
the nominator of $H^{\prime}(x)$ as $J(x)$.

The function $J(x)$ can be rewritten as

$J(x)=x\tilde{x}\left({\sqrt{1+\frac{2}{x}}{{\left({x+1}\right)}^{2}}-\sqrt{1+\frac{2}{{\tilde{x}}}}{{\left({\tilde{x}+1}\right)}^{2}}}\right).$
(68)

Next, we show that $J(x)\leq 0$ when $\frac{{\sqrt{17}-3}}{4}\leq
x\leq\tilde{x}$, and $J(x)>0$ when $x>\tilde{x}$. In specific, let us define

$\vspace{-0.1cm}U(x)\buildrel\Delta\over{=}\sqrt{1+\frac{2}{x}}{\left({x+1}\right)^{2}}.$
(69)

The first-order derivative of $U(x)$ w.r.t. $x$ is given by

$U^{\prime}(x)=\frac{{x+1}}{{{x^{1.5}}\sqrt{x+2}}}\left({2{x^{2}}+3x-1}\right).$
(70)

It can be readily proved that $U^{\prime}(x)\leq 0$ when
$0<x\leq\frac{{\sqrt{17}-3}}{4}$, and $U^{\prime}(x)>0$ when
$x>\frac{{\sqrt{17}-3}}{4}$.
As a result, when $\frac{{\sqrt{17}-3}}{4}\leq x\leq\tilde{x}$, we have
$U(x)<U(\tilde{x})$, and thus $J(x)<0$. From (67), we know that
$H^{\prime}(x)<0$, which means $H(x)$ is a monotonically decreasing function
of $x$. Hence, $H(x)\geq H(0)=0$ holds for $\frac{{\sqrt{17}-3}}{4}\leq
x\leq\tilde{x}$, which equivalently means $F(x)\geq G(x)$. On the other hand,
when $x>\tilde{x}$, we have $U(x)>U(\tilde{x})$, and thus $J(x)>0$, which
means $H^{\prime}(x)>0$ and $H(x)$ is a monotonically increasing function of
$x$. Hence, $H(x)>H(0)=0$ holds for $x\geq\tilde{x}$, and thus $F(x)>G(x)$.

Based on the above analysis, when $x\geq\frac{{\sqrt{17}-3}}{4}$, $F(x)$ is
always no smaller than $G(x)$, which completes the proof.

##  Appendix E Proof of Theorem 3

The first equation in (44) can be readily verified. Then, we focus on the
second equality in (44). The partial derivative of $W({\breve{\bf x}})$ and
$Y({\breve{\bf x}})$ w.r.t. $x_{i}$ are given by

$\frac{{\partial W({\bf{x}})}}{{\partial{x_{i}}}}=\prod\nolimits_{j\neq
i}{(1+{x_{j}})},\frac{{\partial
Y({\bf{x}})}}{{\partial{x_{i}}}}=\lambda{\tau_{i}}x_{i}^{-1}\prod\nolimits_{j\in\mathcal{K}}{x_{j}^{{\tau_{j}}}},i=1,\cdots,K.$
(71)

By substituting ${\bf{x}}=\breve{\bf x}$, $\lambda$ and ${\tau_{i}},\forall i$
into into the above two functions, we can show that

${\left.{\frac{{\partial
W({\bf{x}})}}{{\partial{x_{i}}}}}\right|_{{\bf{x}}={\breve{\bf
x}}}}={\left.{\frac{{\partial
Y({\bf{x}})}}{{\partial{x_{i}}}}}\right|_{{\bf{x}}={\breve{\bf
x}}}},i=1,\cdots,K.$ (72)

Hence, the second equality in (44) is proved.

Now, we begin to prove (42). Before proceeding, we first introduce the
following lemma.

_Lemma 5_ : For any given $\breve{x}\geq 0$, we have the following inequality:

$\frac{{1+x}}{{{x^{\tau}}}}\geq\frac{{1+\breve{x}}}{{{{\breve{x}}^{\tau}}}},$
(73)

where $\tau$ is given by $\tau=\frac{{\breve{x}}}{{1+\breve{x}}}$, and
equality holds only when $x=\breve{x}$.

_Proof_ : The proof is similar to those in Lemma 3, and thus omitted for
simplicity.

By applying inequality (73) for each $x_{i},i=1,\cdots,K$, we have

$\frac{{1+{x_{i}}}}{{x_{i}^{{\tau_{i}}}}}\geq\frac{{1+{{\breve{x}}_{i}}}}{{x_{i}^{{\tau_{i}}}}},i=1,\cdots,K.$
(74)

Then, by multiplying the above $K$ inequalities, we have

$\prod\nolimits_{i\in\mathcal{K}}{\frac{{1+{x_{i}}}}{{x_{i}^{{\tau_{i}}}}}}\geq\prod\nolimits_{i\in\mathcal{K}}{\frac{{1+{x_{i}}}}{{x_{i}^{{\tau_{i}}}}}}$
(75)

which completes the proof.

## References

  * [1] R. Drath and A. Horch, “Industrie 4.0: Hit or hype?[industry forum],” _IEEE industrial electronics mag._ , vol. 8, no. 2, pp. 56–58, 2014.
  * [2] H. Ren, C. Pan, Y. Deng, M. Elkashlan, and A. Nallanathan, “Joint power and blocklength optimization for URLLC in a factory automation scenario,” _IEEE Trans. Wireless Commun._ , 2019.
  * [3] Z. Pang, M. Luvisotto, and D. Dzung, “Wireless high-performance communications: The challenges and opportunities of a new target,” _IEEE Industrial Electronics Mag._ , vol. 11, no. 3, pp. 20–25, 2017.
  * [4] M. Luvisotto, Z. Pang, and D. Dzung, “Ultra high performance wireless control for critical applications: Challenges and directions,” _IEEE Trans. Industrial Informatics_ , vol. 13, no. 3, pp. 1448–1459, 2016.
  * [5] Z. Ma, M. Xiao, Y. Xiao, Z. Pang, H. V. Poor, and B. Vucetic, “High-reliability and low-latency wireless communication for internet of things: Challenges, fundamentals and enabling technologies,” _IEEE Internet Things J._ , 2019.
  * [6] 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, May 2010.
  * [7] C. She, C. Yang, and T. Q. Quek, “Radio resource management for ultra-reliable and low-latency communications,” _IEEE Commun. Mag._ , vol. 55, no. 6, pp. 72–78, 2017.
  * [8] X. Sun _et al._ , “Short-Packet Downlink Transmission With Non-Orthogonal Multiple Access,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 7, pp. 4550–4564, July 2018.
  * [9] Y. Hu _et al._ , “SWIPT-Enabled Relaying in IoT Networks Operating with Finite Blocklength Codes,” _IEEE J. Sel. Areas in Commun._ , pp. 1–1, 2018\.
  * [10] C. Pan, H. Ren, Y. Deng, M. Elkashlan, and A. Nallanathan, “Joint blocklength and location optimization for URLLC-enabled UAV relay systems,” _IEEE Commun. Lett._ , vol. 23, no. 3, pp. 498–501, March 2019.
  * [11] H. Ren, C. Pan, Y. Deng, M. Elkashlan, and A. Nallanathan, “Resource Allocation for Ultra-Reliable and Low-Latency Communications in 5G Mission-Critical IoT Networks,” in _IEEE ICC_ , 2019.
  * [12] J. Chen, L. Zhang, Y.-C. Liang, X. Kang, and R. Zhang, “Resource allocation for wireless-powered IoT networks with short packet communication,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 2, pp. 1447–1461, 2019.
  * [13] C. She, C. Yang, and T. Q. S. Quek, “Joint Uplink and Downlink Resource Configuration for Ultra-Reliable and Low-Latency Communications,” _IEEE Trans. Commun._ , vol. 66, no. 5, pp. 2266–2280, May 2018.
  * [14] W. R. Ghanem, V. Jamali, Y. Sun, and R. Schober, “Resource allocation for multi-user downlink URLLC-OFDMA systems,” _arXiv preprint arXiv:1901.05825_ , 2019.
  * [15] V. K. Huang, Z. Pang, C.-J. A. Chen, and K. F. Tsang, “New trends in the practical deployment of industrial wireless: From noncritical to critical use cases,” _IEEE Industrial Electronics Mag._ , vol. 12, no. 2, pp. 50–58, 2018\.
  * [16] J. G. Andrews, S. Buzzi, W. Choi, S. V. Hanly, A. Lozano, A. C. Soong, and J. C. Zhang, “What will 5G be?” _IEEE J. Sel. Areas Commun._ , vol. 32, no. 6, pp. 1065–1082, 2014.
  * [17] H. V. Cheng, E. Björnson, and E. G. Larsson, “Optimal pilot and payload power control in single-cell massive MIMO systems,” _IEEE Trans. Signal Process._ , vol. 65, no. 9, pp. 2363–2378, 2017.
  * [18] O. Saatlou, M. O. Ahmad, and M. Swamy, “Joint data and pilot power allocation for massive mu-mimo downlink tdd systems,” _IEEE Trans. Circuits and Systems II: Express Briefs_ , vol. 66, no. 3, pp. 512–516, 2019.
  * [19] S. Lu and Z. Wang, “Joint optimization of power allocation and training duration for uplink multiuser mimo communications,” in _2015 IEEE WCNC_. IEEE, 2015, pp. 322–327.
  * [20] T. L. Marzetta, “Noncooperative cellular wireless with unlimited numbers of base station antennas,” _IEEE Trans. Wireless Communs._ , vol. 9, no. 11, pp. 3590–3600, November 2010.
  * [21] J. Shen, J. Zhang, and K. B. Letaief, “Downlink user capacity of massive MIMO under pilot contamination,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 6, pp. 3183–3193, June 2015.
  * [22] M. Mousaei and B. Smida, “Optimizing pilot overhead for ultra-reliable short-packet transmission,” in _IEEE ICC_ , 2017, pp. 1–5.
  * [23] H. Ren, C. Pan, K. Wang, Y. Deng, M. Elkashlan, and A. Nallanathan, “Achievable data rate for URLLC-enabled UAV systems with 3-D channel model,” _IEEE Wireless Commun. Letter,_ , vol. 8, no. 6, pp. 1587–1590, Dec 2019.
  * [24] Z.-Q. Luo and S. Zhang, “Dynamic spectrum management: Complexity and duality,” _IEEE J. Sel. Topics Signal Process._ , vol. 2, no. 1, pp. 57–73, 2008.
  * [25] H. V. Cheng, E. Björnson, and E. G. Larsson, “Uplink pilot and data power control for single cell massive MIMO systems with MRC,” in _IEEE ISWCS_ , 2015, pp. 396–400.
  * [26] O. N. Yilmaz, Y.-P. E. Wang, N. A. Johansson, N. Brahmi, S. A. Ashraf, and J. Sachs, “Analysis of ultra-reliable and low-latency 5G communication for a factory automation use case,” in _IEEE ICCW_ , 2015, pp. 1190–1195.
  * [27] S. Boyd and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [28] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1,” 2014.
  * [29] E. Matskani, N. D. Sidiropoulos, Z.-q. Luo, and L. Tassiulas, “Convex approximation techniques for joint multiuser downlink beamforming and admission control,” _IEEE Trans. Wireless Commun._ , vol. 7, no. 7, pp. 2682–2693, 2008.
  * [30] C. Pan, H. Zhu, N. J. Gomes, and J. Wang, “Joint user selection and energy minimization for ultra-dense multi-channel C-RAN with incomplete CSI,” _IEEE J. Sel. Areas Commun._ , vol. 35, no. 8, pp. 1809–1824, 2017.
  * [31] S. Boyd, S.-J. Kim, L. Vandenberghe, and A. Hassibi, “A tutorial on geometric programming,” _Optimization and engineering_ , vol. 8, no. 1, p. 67, 2007\.
  * [32] S. Shi, M. Schubert, and H. Boche, “Rate optimization for multiuser MIMO systems with linear processing,” _IEEE Trans. Signal Process._ , vol. 56, no. 8, pp. 4020–4030, 2008.
  * [33] M. Chiang, C. W. Tan, D. P. Palomar, D. O Neill, and D. Julian, “Power control by geometric programming,” _Resource Allocation in Next Generation Wireless Networks_ , vol. 5, pp. 289–313, 2005.
  * [34] E. U. T. R. Access, “Further advancements for E-UTRA physical layer aspects,” _3GPP TR 36.814, Tech. Rep._ , 2010.
  * [35] G. W. R, S. Y. Jamali V, and R. Schober, “Resource allocation for multi-user downlink URLLC-OFDMA systems,” _arXiv preprint arXiv:1901.05825, 2019_ , 2016.
  * [36] H. Q. Ngo, E. G. Larsson, and T. L. Marzetta, “Energy and spectral efficiency of very large multiuser MIMO systems,” _IEEE Trans. Commun._ , vol. 61, no. 4, pp. 1436–1449, 2013.
  * [37] A. M. Tulino, S. Verdú _et al._ , “Random matrix theory and wireless communications,” _Foundations and Trends® in Commun. and Inf. Theory_ , vol. 1, no. 1, pp. 1–182, 2004.