# RIS-assisted Cell-Free MIMO with Dynamic Arrivals and Departures of Users: A
Novel Network Stability Approach

Charbel Bou Chaaya1, Mohamad Assaad2 and Tijani Chahed1
1 Institut Polytechnique de Paris, Telecom SudParis, 19 Place Marguerite
Perey, 91120 Palaiseau, France
2 University of Paris-Saclay, CNRS, CentraleSupelec, Laboratoire des signaux
et systèmes, 91190 Gif-sur-Yvette, France

###### Abstract

Reconfigurable Intelligent Surfaces (RIS) have recently emerged as a hot
research topic, being widely advocated as a candidate technology for next
generation wireless communications. These surfaces passively alter the
behavior of propagation environments enhancing the performance of wireless
communication systems. In this paper, we study the use of RIS in cell-free
multiple-input multiple-output (MIMO) setting where distributed service
antennas, called Access Points (APs), simultaneously serve the users in the
network. While most existing works focus on the physical layer improvements
RIS carry, less attention has been paid to the impact of dynamic arrivals and
departures of the users. In such a case, ensuring the stability of the network
is the main goal. For that, we propose an optimization framework of the phase
shifts, for which we derived a low-complexity solution. We then provide a
theoretical analysis of the network stability and show that our framework
stabilizes the network whenever it is possible. We also prove that a low
complexity solution of our framework stabilizes a guaranteed fraction (higher
than 78.5%) of the stability region. We provide also numerical results that
corroborate the theoretical claims.

###### Index Terms:

Reconfigurable Intelligent Surfaces (RIS), Cell-free MIMO, Stability Region,
Flow-level modeling, Dynamic user configuration

##  I INTRODUCTION

Unprecedented communication requirements for Beyond 5G networks will call for
essentially new schemes. Numerous technologies have been investigated in the
last few decades to satisfy the exponential surge for wireless connectivity,
notably massive Multiple-Input Multiple-Output (mMIMO)[1, 2], ultra dense
networks (UDN) and millimeter waves (mmWave) [3]. The reliable service offered
by these innovations however, comes at the cost of expensive hardware and high
energy consumption. Consequently, extensive research is ongoing to find novel
designs for sustainable wireless networks with spectrum and power
efficiencies, and low hardware cost.

Reconfigurable intelligent surfaces (RIS), or intelligent reconfigurable
surfaces (IRS), have recently emerged as a paradigm that can leverage
engineered scattering surfaces to transmit and receive information [4]. As a
candidate for 6G communications, these surfaces are capable of shaping the
radio waves at the electromagnetic level, without the need of radio frequency
chains, power amplifiers, or signal processing in general. Specifically, a RIS
is a planar array composed of a large number of passive scattering elements
(such as printed dipoles), that induce controlled phase shifts on the incident
signal. These phase shifts are independent from the impinging wave, and thus
the RIS elements collaborate to modify the propagation of the reflected
signal. Eventually, by smartly tuning the phase shift of each element, the
reflected signals from different paths can be coherently combined at the
desired receiver to improve the quality of the received signal. Accordingly,
the RIS can maneuver the propagation environment to boost coverage over the
area of interest and avoid blockages.

Although incorporating RIS in wireless systems entails challenging
optimization problems, its usage has shown real potential in enhancing point-
to-point communications, as well as uplink and downlink multi-user schemes
[5]. Within this framework, recent works on RIS-aided wireless systems have
predominantly studied the physical layer advantages, while the protocol
aspects for their integration have been generally overlooked. In addition,
less attention has been paid to control signaling and network stability,
defined as the set of users arrival rates that the network can serve within
finite time, which are imperative for their operation. For instance, the
authors in [6] conceived a simple random access algorithm for RIS-assisted
communications, and [7] presented a next generation multiple access (NGMA)
scheme with RIS. Moreover, [8] tackled the problem of non-orthogonal multiple
access (NOMA) with RIS, and [9] proposed a resource allocation scheme with
rate splitting multiple access (RSMA). On the other hand, [10] studied a RIS-
aided cell-free MIMO system where distributed access points (APs) in the
coverage area cooperate via a backhaul unit to serve the users. This
architecture alleviates the intercell interference that characterizes mMIMO
[11], and will be used as the physical layer in this work. Furthermore, [12]
considers a RIS-assisted cell free massive MIMO scheme with random phase
shifts at the RIS and conjugate beamforming at the massive base station (BS)
antennas. A closed-form expression of the achievable rate is obtained in [12]
and a gain in terms of system coverage and user rate is shown. The
aforementioned works focused on the physical layer without considering the
impact of dynamic traffic.

Contrary to the previous studies that consider a fixed set of users in the
system, we focus in this paper on a more realistic scenario where users
dynamically connect to the network, exchange bursts of information and then
leave. We first rigorously describe the physical layer of the communication
and the traffic/flow level, where flows represent file transfers. The primary
difficulty of this problem is twofold: firstly the presence of the RIS that
involves complex analytical forms for the channel modelling, and secondly the
dynamic nature of the traffic and the interference between a variable number
of users.

In the context of dynamic arrivals departures of users, the goal is to ensure
that the network is stable [13]. Network stability implies that all
files/users will be served in a finite time. Network stability with dynamic
population has been considered in [13] to deal with power control in massive
MIMO. In this paper, we consider a different problem of phase shifts in RIS-
assisted cell-free MIMO. Instead of considering a sum rate maximization, we
show that achieving network stability requires to solve a phase shift
framework that depends on the current flows intensity and
$\operatorname{SINR}$ (Signal to Interference and Noise Ratio) of the users.
Typically, globally optimal phase shift solutions cannot be obtained for RIS-
related problems due to their non-convexity. We show that a low complexity
solution of our proposed framework achieves a guaranteed fraction $\geq
78.5\%$ of the network stability region.

To the best of our knowledge, this is the first work that considers a RIS-
assisted wireless environment with dynamic arrival and departures of the
users. The remainder of this paper is organized as follows. In section II, the
system model is provided with the necessary definitions of stability. The
targeted phase shift design problem is then inspected in section III, and the
stability analysis is presented in section IV. Numerically evaluated results
are given in section V, while section VI concludes the paper.

##  II SYSTEM MODEL

###  II-A CHANNEL MODEL

We consider the uplink of a RIS-aided cell-free MIMO system illustrated in
Figure 1, where $N$ access points (APs) equipped with a single antenna each,
serve $K$ single antenna users who transmit their signals simultaneously on
the same time-frequency resource. All APs are linked to a central processing
unit (CPU) via a backhaul network. The communication is assisted by an RIS
comprising of $M$ reflecting elements. Typically, the number of RIS engineered
elements $M$ is very large. We focus on the scenario where the direct channels
between the users and the APs are blocked, and the communication is ensured
only by the RIS, that can modify the phases of incident signals. These phase
shifts introduced by the RIS are adjusted adaptively by a controller linked to
the CPU.

We assume a block fading model where the channel realizations are generated
randomly and are independent between blocks. The channels between user
$k=1,\dots,K$ and the RIS, and the RIS and AP $n=1,\dots,N$ are respectively
denoted by $\bm{g}_{k}\in\mathbb{C}^{M\times 1}$ and
$\bm{h}_{n}\in\mathbb{C}^{M\times 1}$, and are modeled as ergodic processes
with independent realizations $\bm{g}_{k}\sim\mathcal{CN}\left(\bm{0}_{M\times
1},\,\overline{\bm{R}}_{k}\right)$ and
$\bm{h}_{n}\sim\mathcal{CN}\left(\bm{0}_{M\times
1},\,\widetilde{\bm{R}}_{n}\right)$, where
$\overline{\bm{R}}_{k},\widetilde{\bm{R}}_{n}\in\mathbb{C}^{M\times M}$ are
deterministic positive semi-definite covariance matrices that describe the
spatial correlation between the channels of the RIS elements. Further, we can
write $\overline{\bm{R}}_{k}=\overline{\alpha}_{k}\bm{R}_{\text{r}}$ and
$\widetilde{\bm{R}}_{n}=\widetilde{\alpha}_{n}\bm{R}_{\text{t}}$, with
$\overline{\alpha}_{k}$, $\widetilde{\alpha}_{n}\in\mathbb{C}$ being the the
large-scale fading coefficients, and $\bm{R}_{\text{r}}$, $\bm{R}_{\text{t}}$
being the RIS receive and transmit correlation matrices respectively. These
matrices can follow a general model that depends on the size, distance and
layout of RIS scattering elements, and the propagation environment. For
example, a popular model for such planar arrays is the Kroenecker model where
the correlation matrix is Hermitian Toeplitz with exponential entries
$\rho^{\absolutevalue{j-i}}$ at row $i$ and column $j$, if $j\geq i$, where
$\rho\in\mathbb{C}$ is the correlation coefficient that satisfies
$\absolutevalue{\rho}\leq 1$. For the case of isotropic scattering in front of
the RIS, the correlation matrix can be explicitly obtained from [14,
Proposition 1]. We assume for the rest of this paper, that the trasmit and
receive correlation matrices are not equal, and the only technical assumption
is that their diagonal elements are equal to unity due to power normalization.

The cascaded channel between user $k$ and AP $n$ is then written:

$u_{n,k}=\bm{h}_{n}^{\text{H}}\,\bm{\Theta}\,\bm{g}_{k}$ (1)

where $\bm{\Theta}=\operatorname{diag}(\bm{\phi})$ is the RIS reflection
matrix, with $\bm{\phi}=\left[\phi_{1},\dots,\phi_{M}\right]^{\intercal}$,
$\phi_{m}=e^{j\theta_{m}}$ and
$\bm{\theta}=\left[\theta_{1},\dots,\theta_{M}\right]^{\intercal}\in\left[0,2\pi\right]^{M}$
are the phase shifts of the $M$ elements. We define the aggregated channel
from user $k$ to the APs as $\bm{u}_{k}=[u_{1,k},\dots,u_{N,k}]^{\intercal}$,
where $u_{n,k}=\sum_{m=1}^{M}h_{nm}\,g_{km}\,e^{j\theta_{m}}$ is the virtual
link between the $k^{\text{th}}$ user and the $n^{\text{th}}$ AP.

Figure 1: System Model

Similar models have been studied for massive MIMO systems in [15], and for
cell-free MIMO in [10]. For instance, the second moment of the aggregated
channel can be directly obtained as:

$\displaystyle\begin{split}\mathbb{E}\left[\absolutevalue{u_{n,k}}^{2}\right]&=\mathbb{E}\left[\absolutevalue{\bm{h}^{\text{H}}_{n}\,\bm{\Theta}\,\bm{g}_{k}}^{2}\right]\\\
&=\mathbb{E}\left[\tr\left(\bm{h}_{n}^{\text{H}}\,\bm{\Theta}\,\bm{g}_{k}\,\bm{g}^{\text{H}}_{k}\,\bm{\Theta}^{\text{H}}\,\bm{h}_{n}\right)\right]\\\
&=\tr\left(\bm{\Theta}^{\text{H}}\,\mathbb{E}\left[\bm{h}_{n}\,\bm{h}_{n}^{\text{H}}\right]\,\bm{\Theta}\,\mathbb{E}\left[\bm{g}_{k}\,\bm{g}^{\text{H}}_{k}\right]\right)\\\
&=\tr\left(\bm{\Theta}^{\text{H}}\,\widetilde{\bm{R}}_{n}\,\bm{\Theta}\,\overline{\bm{R}}_{k}\right).\end{split}$
(2)

We can also see that the channels between two different users and two
different APs are mutually independent, $\forall n\neq n^{\prime},k\neq
k^{\prime},\mathbb{E}\left[u_{n,k}\,u_{n^{\prime},k^{\prime}}^{*}\right]=\mathbb{E}\left[u_{n,k}\right]\,\mathbb{E}\left[u_{n^{\prime},k^{\prime}}^{*}\right]=0$.
Moreover, the channels between a user and two different APs are uncorrelated,
$\forall n\neq
n^{\prime},\mathbb{E}\left[u_{n,k}\,u_{n^{\prime},k}^{*}\right]=\mathbb{E}\left[\left(\bm{h}_{n}^{\text{H}}\bm{\Theta}\bm{g}_{k}\right)\left(\bm{h}_{n^{\prime}}^{\text{H}}\bm{\Theta}\bm{g}_{k}\right)^{*}\right]=0$,
and the channels between two different users and the same AP are also
uncorrelated, $\forall k\neq
k^{\prime},\mathbb{E}\left[u_{n,k}\,u_{n,k^{\prime}}^{*}\right]=\mathbb{E}\left[\left(\bm{h}_{n}^{\text{H}}\bm{\Theta}\bm{g}_{k}\right)\left(\bm{h}^{\text{H}}_{n}\bm{\Theta}\bm{g}_{k^{\prime}}\right)^{*}\right]=0$.
Finally, we define the aggregated channel covariance matrix
$\bm{Q}_{k}=\mathbb{E}\left[\bm{u}_{k}\,\bm{u}_{k}^{\text{H}}\right]=\overline{\alpha}_{k}\tr\left(\bm{\Theta}^{\text{H}}\,\bm{R}_{\text{t}}\,\bm{\Theta}\,\bm{R}_{\text{r}}\right)\operatorname{diag}\left(\widetilde{\alpha}_{1},\dots,\widetilde{\alpha}_{N}\right)$.

###  II-B UPLINK TRANSMISSION

In a simultaneous manner, the $K$ users transmit their signals to the $N$ APs,
that know an estimate $\widehat{u}_{n,k}$ of their channels, computed during
the training phase. The baseband received signal at AP $n$ is:

$y_{n}=\sum_{k=1}^{K}u_{n,k}\,s_{k}+n_{n}$ (3)

where the additive noise at the $n^{\text{th}}$ AP
$n_{n}\sim\mathcal{CN}\left(0,\,\sigma^{2}\right)$ and the uplink transmitted
symbol $s_{k}\in\mathbb{C}$ satisfies
$\mathbb{E}\left[\absolutevalue{s_{k}}^{2}\right]=P_{k}$, with $P_{k}>0$ being
the transmit power budget of user $k$.

For data detection of the $k^{\text{th}}$ symbol, the $n^{\text{th}}$ AP
multiplies the received signal $y_{n}$ with the conjugate of its locally
obtained channel estimate. The acquired terms are then sent to the CPU via the
backhaul network. Therefore, the decision statistic for user $k$’s symbol
$s_{k}$ reads:

$r_{k}=\sum_{n=1}^{N}\widehat{u}^{*}_{n,k}\,y_{n}=\sum_{n=1}^{N}\sum_{j=1}^{K}\widehat{u}^{*}_{n,k}\,u_{n,j}s_{j}+\sum_{n=1}^{N}\widehat{u}^{*}_{n,k}\,n_{n}.$
(4)

In the following, to maintain analytical tractability, we assume that all
users employ mutually orthogonal reverse link pilot sequences. Accordingly,
channel estimation error can be neglected when the number of RIS elements $M$
is large, as shown in [16]. Hence, $\widehat{u}_{n,k}=u_{n,k}$.

By utilizing the use-and-then-forget bounding technique [17, Theorem 4.4] and
the received signal expression in (4), we can show that an achievable uplink
rate expression for user $k$ is $\log\left(1+\operatorname{SINR}_{k}\right)$,
where the effective uplink $\operatorname{SINR}$ is:

$\operatorname{SINR}_{k}=\frac{\absolutevalue{\mathtt{DS}_{k}}^{2}P_{k}}{\sum\limits_{\begin{subarray}{c}j=1\\\
j\neq
k\end{subarray}}^{K}\mathbb{E}\left[\absolutevalue{\mathtt{UI}_{j,k}}^{2}\right]P_{j}+\mathbb{E}\left[\absolutevalue{\mathtt{BU}_{k}}^{2}\right]P_{k}+\mathbb{E}\left[\absolutevalue{\mathtt{NO}_{k}}^{2}\right]}$
(5)

where we have

$\absolutevalue{\mathtt{DS}_{k}}^{2}=\absolutevalue{\mathbb{E}\left[\sum\limits_{n=1}^{N}\widehat{u}^{*}_{n,k}u_{n,k}\right]}^{2}=\absolutevalue{\mathbb{E}\left[\bm{u}_{k}^{\text{H}}\bm{u}_{k}\right]}^{2}=\left(\tr\left(\bm{Q}_{k}\right)\right)^{2}$
(6)
$\displaystyle\begin{split}\mathbb{E}\left[\absolutevalue{\mathtt{UI}_{j,k}}^{2}\right]&=\mathbb{E}\left[\absolutevalue{\sum\limits_{n=1}^{N}\widehat{u}^{*}_{n,k}\,u_{n,j}}^{2}\right]=\mathbb{E}\left[\absolutevalue{\bm{u}_{k}^{\text{H}}\bm{u}_{j}}^{2}\right]\\\
&=\tr\left(\bm{Q}_{k}\,\bm{Q}_{j}\right)\end{split}$ (7)
$\mathbb{E}\left[\absolutevalue{\mathtt{NO}_{k}}^{2}\right]=\mathbb{E}\left[\absolutevalue{\sum\limits_{n=1}^{N}\widehat{u}^{*}_{n,k}\,n_{n}}^{2}\right]=\mathbb{E}\left[\absolutevalue{\bm{u}_{k}^{\text{H}}\bm{n}}^{2}\right]=\sigma^{2}\tr\left(\bm{Q}_{k}\right)$
(8)
$\displaystyle\begin{split}\mathbb{E}\left[\absolutevalue{\mathtt{BU}_{k}}^{2}\right]&=\mathbb{E}\left[\absolutevalue{\sum\limits_{n=1}^{N}\widehat{u}^{*}_{n,k}\,u_{n,k}-\mathbb{E}\left[\sum\limits_{n=1}^{N}\widehat{u}^{*}_{n,k}\,u_{n,k}\right]}^{2}\right]\\\
&=\mathbb{V}\left[\bm{u}_{k}^{\text{H}}\bm{u}_{k}\right]=\tr\left(\bm{Q}_{k}^{2}\right)\end{split}$
(9)

where the final equality follows from [17, Lemma B.14], because the aggregated
channels can be approximated by Gaussian distributions when $M$ is large as
shown in [16]. These variables represent the strength of the desired signal
$\left(\mathtt{DS}_{k}\right)$, the beamforming gain uncertainty
$\left(\mathtt{BU}_{k}\right)$, the interference caused by user $j$ on user
$k$ $\left(\mathtt{UI}_{j,k}\right)$ and the power of additive noise
$\left(\mathtt{NO}_{k}\right)$. After some algebraic manipulations, we
simplify the $\operatorname{SINR}$ expression to:

$\operatorname{SINR}_{k}=\frac{\overline{\alpha}_{k}^{2}\left(\sum_{n=1}^{N}\widetilde{\alpha}_{n}\right)^{2}\eta\,P_{k}}{\sum_{j=1}^{K}\overline{\alpha}_{k}\,\overline{\alpha}_{j}\left(\sum_{n=1}^{N}\widetilde{\alpha}_{n}^{2}\right)\eta\,P_{j}+\sigma^{2}\,\overline{\alpha}_{k}\left(\sum_{n=1}^{N}\widetilde{\alpha}_{n}\right)}$
(10)

where
$\eta=\tr\left(\bm{\Theta}^{\text{H}}\,\bm{R}_{\text{t}}\,\bm{\Theta}\,\bm{R}_{\text{r}}\right)$.

Figure 2: Flows Arrival Model

###  II-C USERS / FLOWS ARRIVAL MODEL

All previous works on RIS-aided systems make the assumption that the users are
static in the network, in the sense that a fixed number of devices constantly
communicate with the BS or the APs. Nevertheless, in real settings, users join
the network dynamically, exchange bursts of data with the APs, and leave once
they are served. Subsequently, we adopt the model used in [13] for power
control in massive MIMO networks, and consider this dynamic user population
around the RIS.

Specifically, we consider a finite, but possibly large, number of locations in
the network in the vicinity of the RIS. Accordingly, we use $K$ to denote the
number of locations instead of users. Let $X_{k}(t)$ represent the number of
users at location $k$ at time $t$. Each user has one flow to be served, where
a flow depicts a file that the user would like to convey to the APs. Thus, the
words ‘flow’ and ‘user’ indicate a certain file transfer, and will be
interchanged for the rest of this paper. Note that a similar model can be used
for the downlink. Also, one can notice that a model where each user has
multiple flows is a simple extension to this case. The arrival of the flows at
location $k$ is modeled as a Poisson process with rate $\lambda_{k}$. These
rates represent the average number of users arriving to each location. The
size of the files to be transmitted is an exponentially distributed random
variable $S_{k}$ at location $k$ with mean $\mathbb{E}\left[S_{k}\right]$. The
$\operatorname{SINR}$ serving the flow at location $k$ is obtained from (10).
Figure 2 displays the considered flows arrival model.

We consider distinct time-scales for the physical layer and the flow level,
and we use the time index $t$ to refer to the time at the level of flows. In
practice, multiple physical layer timeslots occur between times $t$ and $t+1$.
In the rest of this paper, we will suppress the use of the time index when
ambiguity is unlikely.

Next, we define the network stability to be used in the stability analysis.

######  Definition 1 (Strong Stability [18]).

The network is said to be stable if
$\lim\limits_{T\to\infty}\sup\,\frac{1}{T}\sum_{T=0}^{T-1}\sum_{k=1}^{K}\mathbb{E}\left[X_{k}(t)\right]<\infty$

In simple terms, based on this stability definition, if the network is stable
then all arriving users are expected to transmit their files in finite time.
We let $\bm{\lambda}=\left[\lambda_{1},\dots,\lambda_{K}\right]^{\intercal}$
denote the vector of user arrival rates. We can now present the definition of
a stability region.

######  Definition 2.

The stability region is the set of all mean arrival rate vectors for which
there exists a RIS configuration that makes the network stable.

We are now interested in finding a RIS configuration that stabilizes the
network whenever possible.

##  III RIS PHASE SHIFTS DESIGN

First, we formulate the problem we are interested in solving. Let
$\bm{\Lambda}^{\text{max}}$ be the stability region under RIS configurations.
Our target is then:

Choose $\displaystyle\bm{\theta}\in\left[0,2\pi\right]^{M}$ such that
$\displaystyle\forall\bm{\lambda}\in\bm{\Lambda}^{\textrm{max}},\,\textrm{the
network is stable.}$

Note that the flows’ arrival rates $\lambda_{k}$ might not be known a priori
at the CPU; and hence, the network must be stabilized only by modifying the
RIS phase shifts for any arrival rates lying inside the stability region. To
accomplish this, we consider the following optimization problem:

(P1) $\displaystyle\underset{\bm{\theta}}{\textrm{maximize}}\quad$
$\displaystyle
f\left(\bm{\theta}\right)=\sum_{k=1}^{K}X_{k}\,U_{k}\left(\bm{\theta}\right)$
subject to $\displaystyle\quad 0\leq\theta_{i}\leq 2\pi\quad i=1,\dots,M$

where we choose the utility functions:
$U_{k}\left(\bm{\theta}\right)=\log\left(\operatorname{SINR}_{k}\right)$111Throughout
this paper, we use $\log\left(\cdot\right)$ as the natural logarithm function
for simplicity, even in the bitrate expression. The analysis remains true for
any other base, since all logarithms are equal up to a scaling factor.. This
choice is motivated by the concavity of the equivalent utility functions with
respect to the allocated rates [13]. This asset will be discussed in the next
section.

######  Remark 1.

We strongly emphasize that we do not make any kind of approximations in (P1),
such as
$\log\left(1+\operatorname{SINR}\right)\approx\log\left(\operatorname{SINR}\right)$
for high $\operatorname{SINR}$ values. The objective function can be seen as a
sum of each location’s utility functions, selected as
$\log\left(\operatorname{SINR}_{k}\right)$, weighted by the traffic volume at
each location $X_{k}$. Regardless, the rate is still equal to
$\log\left(1+\operatorname{SINR}\right)$. Informally, (P1) can be seen as a
variation to the proportional fairness problem but at the
$\operatorname{SINR}$ level.

Solving this optimization problem is challenging, because of the intricate
form of its objective function. To mitigate this difficulty, we start by
showing that the objective function $f\left(\bm{\theta}\right)$ is increasing
in $\eta$:

$\displaystyle\begin{split}&\frac{\partial
f\left(\bm{\theta}\right)}{\partial\,\eta}=\sum_{k=1}^{K}X_{k}\times\\\
&\frac{\sigma^{2}\,\overline{\alpha}_{k}\left(\sum_{n=1}^{N}\widetilde{\alpha}_{n}\right)}{\eta\left(\sigma^{2}\,\overline{\alpha}_{k}\left(\sum_{n=1}^{N}\widetilde{\alpha}_{n}\right)+\sum_{j=1}^{K}\overline{\alpha}_{k}\,\overline{\alpha}_{j}\left(\sum_{n=1}^{N}\widetilde{\alpha}_{n}^{2}\right)\eta\,P_{j}\right)}>0.\end{split}$
(11)

On the other hand, notice that222This is a direct application of [19, Theorem
1.11]. $\eta=\tr\left(\bm{\Theta}^{\text{H}}\,\bm{R}_{\text{
t}}\,\bm{\Theta}\,\bm{R}_{\text{r}}\right)=\bm{\phi}^{\text{H}}\,\bm{R}\,\bm{\phi}$,
where $\bm{R}=\bm{R}_{\text{ t}}\odot\bm{R}^{\intercal}_{\text{r}}$ is
positive semi-definite by virtue of the Schur product theorem [20]. Thus, (P1)
is equivalent to:

(P2) $\displaystyle\underset{\bm{\phi}}{\textrm{maximize}}\quad$
$\displaystyle\bm{\phi}^{\text{H}}\,\bm{R}\,\bm{\phi}$ subject to
$\displaystyle\absolutevalue{\phi_{i}}=1$ $\displaystyle\quad i=1,\dots,M$

(P2) is a complex quadratic optimization problem, that is in the class of NP-
hard problems; and therefore a globally optimal solution cannot be obtained in
polynomial time. However, we know in this case, that applying the semi-
definite relaxation (SDR) technique reported in [21] can bound the error
committed while approximating the solution. As such, the following semi-
definite program (SDP) provides a relaxation for (P2):

(P3) $\displaystyle\underset{\bm{\Phi}}{\textrm{maximize}}\quad$
$\displaystyle\tr\left(\bm{R}\,\bm{\Phi}\right)$ subject to
$\displaystyle\bm{\Phi}\succeq\bm{0},\;\bm{\Phi}_{i,i}=1$ $\displaystyle\quad
i=1,\dots,M$

where $\bm{\Phi}=\bm{\phi}\,\bm{\phi}^{\text{H}}$. (P3) is standard convex
optimization problem that can be optimally solved by invoking interior-point
based solvers, such as CVX [22]. If the obtained solution $\bm{\Phi}^{\star}$
is rank-one, then it is also the optimal solution to (P2). But since this is
not generally the case, we construct a rank-one solution to (P2) from
$\bm{\Phi}^{\star}$ using Guassian randomization as follows. We first compute
the eigenvalue decomposition
$\bm{\Phi}^{\star}=\bm{U}\,\bm{\Sigma}\,\bm{U}^{\text{H}}$, where $\bm{U}$ and
$\bm{\Sigma}$ are a unitary and a diagonal matrix respectively. Then we obtain
a sub-optimal solution to (P2) as
$\overline{\bm{\phi}}=\bm{U}\,\bm{\Sigma}^{1/2}\,\bm{r}$ where
$\bm{r}\in\mathbb{C}^{M\times 1}$ is a random vector generated according to
$\bm{r}\sim\mathcal{CN}\left(\bm{0}_{M\times 1},\,\bm{I}_{M}\right)$. By
drawing independent random vectors $\bm{r}$, the objective function of (P2) is
approximated by the maximum one attained by the corresponding
$\overline{\bm{\phi}}$. We finally recover a suboptimal solution to (P2) as a
vector $\bm{\phi}$ with elements
$\phi_{i}=e^{j\arg\left(\overline{\phi}_{i}\right)},1\leq i\leq M$. It has
been shown in [21], [23] that such techniques with a sufficiently large number
of randomizations, guarantee a $\frac{\pi}{4}$-approximation to (P2), in
following sense:

$\frac{\pi}{4}\,\eta^{\text{opt}}\leq\eta^{\text{sub}}\leq\eta^{\text{opt}}\quad\text{or}\quad\frac{\pi}{4}\leq\gamma=\frac{\eta^{\text{sub}}}{\eta^{\text{opt}}}\leq
1$ (12)

where $\eta^{\text{opt}}$ and $\eta^{\text{sub}}$ are the optimal and
approximate solutions of (P2), and $\gamma$ is the approximation accuracy.

It is worth noting that a similar method can be applied to solve the same
problem in the case the RIS phase shifts can only take a finite number of
discrete values. In such scenarios,
$\bm{\phi}\in\mathcal{V}_{L}=\\{1,\omega,\dots,\omega^{L-1}\\}$, where
$\omega$ is the principal $L^{\text{th}}$ root of unity and $L$ is the number
of phase shift levels. The minimum approximation accuracy becomes
$\frac{\left(L\sin\left(\frac{\pi}{L}\right)\right)^{2}}{4\pi}$.

Due to the complexity of the channel induced by the presence of a RIS in the
network, (P1) cannot be optimally solved. This will affect the stability
analysis which we will discuss later. Yet, we could at least bound the error
carried out with the aforementioned solution. This error will however decrease
the stability region insured by solving (P1) to a smaller region. Nonetheless,
we can still acquire the closed-form solution to (P1) in the special case
where the transmit and receive correlation matrices are equal.

######  Corollary 1.

If the transmit and receive correlation matrices are equal
$\bm{R}_{\text{t}}=\bm{R}_{\text{r}}$, any equal RIS phase shift configuration
$\left(\theta_{1}=\theta_{2}=\dots=\theta_{M}\right)$ is the optimal maximizer
of (P1).

###### Proof.

Notice that since $\bm{R}_{\text{t}}$ and $\bm{R}_{\text{r}}$ are Hermitian,
$\bm{R}=\bm{R}_{\text{t}}\odot\bm{R}^{\intercal}_{\text{r}}$ is symmetric, and
has real non negative elements. Then,

$\displaystyle\begin{split}\eta&=\tr\left(\bm{\Theta}^{\text{H}}\,\bm{R}_{\text{t}}\,\bm{\Theta}\,\bm{R}_{\text{r}}\right)=\sum_{m=1}^{M}\sum_{n=1}^{M}\left[\bm{R}\right]_{m,n}e^{j\theta_{n}-j\theta_{m}}\\\
&=\sum_{m=1}^{M}\sum_{n=1}^{M}\left[\bm{R}\right]_{m,n}\cos\left(\theta_{n}-\theta_{m}\right)\end{split}$
(13)

which is maximized for any equal phase shifts configuration. ∎

When the correlation matrices are equal, Corollary 1 provides a simple yet
effective option for the phase shifts design, and optimally solves the desired
problem.

##  IV STABILITY ANALTSIS OF THE NETWORK

In this part, we will show that the RIS configuration presented in the
previous section, i.e. the one that solves (P1), guarantees the network
stability for any flow arrival rates inside the stability region
$\bm{\Lambda}^{\text{max}}$ when (P1) is solved with no optimality gap, and
for all flow arrival rates that are within a factor proportional to the
optimality gap otherwise. To do so, we utilize fluid limit analysis to
demonstrate the stability. In other terms, for the stochastic process
$\mathbf{X}(t)=\left[X_{1}(t),\dots,X_{K}(t)\right]^{\intercal}$ that portrays
the flows’ volumes, we introduce a deterministic process $\mathbf{Y}(t)$ that
approximates the evolution of $\mathbf{X}$ subject to a certain limiting
criteria. It is known that if the fluid limit associated to a stochastic
process reaches zero in finite time, then the process itself is stable.
Furthermore, to show that $\mathbf{Y}$ is stable, it is sufficient to prove
that its corresponding Lyapunov function exhibits negative drift under the
selected RIS configuration [24], [25].

Recall that the users’ arrivals follow a Poisson distribution with mean
$\lambda_{k}$ at location $k$. We suppose that when a new user joins the
network, he directly starts a connection with the APs and transmits a file
having an average size of $\mathbb{E}\left[S_{k}\right]=1\,\forall\,k$. It
will then be straightforward to extend this situation to different file size
means and renewal arrival processes. In this context, it is clear that
$\mathbf{X}$ is a Markov process that has the following evolution at each
timeslot $t$, and at each location $k$:

$\displaystyle X_{k}\longrightarrow X_{k}+1$ at rate
$\displaystyle\lambda_{k}$ $\displaystyle X_{k}\longrightarrow X_{k}-1$ at
rate $\displaystyle R_{k}$

where $R_{k}$ is the physical layer rate allocated to location $k$ between
time instants $t$ and $t+1$. This rate is expressed as
$R_{k}=\log\left(1+\operatorname{SINR}_{k}\right)$, where the
$\operatorname{SINR}$ at each location is given by (10).

After introducing all the necessary ingredients, we now state the main result
of this paper.

######  Theorem 1.

By choosing the RIS phase shifts that solve (P1), the network is stable for
any user arrival rate vector $\gamma\bm{\lambda}$, where
$\bm{\lambda}\in\bm{\Lambda}^{\text{max}}$ and
$\gamma\in\left[\frac{\pi}{4},1\right]$ is the approximation accuracy of (P2).

The aforementioned theorem implies that our proposed solution allows achieving
a guaranteed fraction higher than $\pi/4$ (i.e. higher than 78.5%) of the
stability region.

###### Proof.

Before proving the theorem, we provide a useful lemma for the demonstration.

######  Lemma 1.

$\forall\,\gamma\in(0,1],\forall\,x>0$, $\log\left(1+\gamma
x\right)\geq\gamma\log\left(1+x\right)$, with equality when $\gamma=1$.

###### Proof.

When $\gamma=1$ the equality is obvious. Otherwise,
$\forall\,\gamma\in(0,1),\forall\,x>0$, the proof follows by composing both
sides of Bernoulli’s inequality [26] $\left(1+x\right)^{\gamma}<1+\gamma x$
with the increasing function $x\mapsto\log(x)$. ∎

The proof consists of studying the fluid system obtained when the initial
number of flows grows to infinity. Specifically, we consider the set of fluid
limits defined by:

$Y_{k}(t)=\lim_{\beta\to\infty}\frac{X_{k}\left(\beta
t\right)}{\beta}\quad\text{with}\quad\sum_{k=1}^{K}X_{k}(0)=\beta.$

Notice that if the limit exists, $\sum\limits_{k=1}^{K}Y_{k}(0)=1$. Given this
initial distribution of the fluid system $\mathbf{Y}(0)$, the evolution of
$\mathbf{Y}(t)$ is uniquely defined by the following set of differential
equations, given by the strong law of large numbers:

$\derivative{t}Y_{k}(t)=\lambda_{k}-R_{k}\quad\text{for all}\;k,t\;\text{such
that}\,Y_{k}(t)>0.$ (14)

Now, we define $\mathbf{R}=\left[R_{1},\dots,R_{K}\right]^{\intercal}$ as a
vector containing the allocated rates at each location. Recall that the rate
is calculated by $R_{k}=\log\left(1+\operatorname{SINR}_{k}\right)$.
Conversely we can write
$\log\left(\operatorname{SINR}_{k}\right)=\log\left(e^{R_{k}}-1\right)$. Then,
(P1) is equivalent to the following optimization problem:

(P4) $\displaystyle\underset{\mathbf{R}}{\textrm{maximize}}\quad$
$\displaystyle\sum_{k=1}^{K}Y_{k}\log\left(e^{R_{k}}-1\right)$ subject to
$\displaystyle 0\leq f_{i}\left(\mathbf{R}\right)\leq 2\pi,\quad i=1,\dots,M$

where $\mathbf{X}$ is interchanged with its limit $\mathbf{Y}$, and the
constraint functions $f_{i}\left(\mathbf{R}\right)$ are the equivalent of the
phase shift constraints on $\bm{\theta}$. Let
$\operatorname{SINR}_{k}\left(\eta\right)$ denote the $\operatorname{SINR}$
value at location $k$ for a given value of $\eta$. Since
$\operatorname{SINR}_{k}\left(\eta\right)$ increases with $\eta$, we have from
(12):

$\displaystyle\begin{split}\gamma\,\operatorname{SINR}_{k}\left(\eta^{\text{opt}}\right)\overset{\text{(a)}}{\leq}\operatorname{SINR}_{k}\left(\gamma\,\eta^{\text{opt}}\right)&\leq\operatorname{SINR}_{k}\left(\eta^{\text{sub}}\right)\\\
&\leq\operatorname{SINR}_{k}\left(\eta^{\text{opt}}\right)\end{split}$ (15)

where (a) is a straightforward bound due to $0<\gamma\leq 1$. Moreover, for
any $x>0$ the function $x\mapsto\log\left(1+x\right)$ is increasing and thus,

$\displaystyle\begin{split}\gamma\,\log\left(1+\operatorname{SINR}_{k}\left(\eta^{\text{opt}}\right)\right)\overset{\text{(b)}}{\leq}\log\left(1+\gamma\,\operatorname{SINR}_{k}\left(\eta^{\text{opt}}\right)\right)\\\
\leq\underbrace{\log\left(1+\operatorname{SINR}_{k}\left(\eta^{\text{sub}}\right)\right)}_{R_{k}^{\text{sub}}}\leq\underbrace{\log\left(1+\operatorname{SINR}_{k}\left(\eta^{\text{opt}}\right)\right)}_{R_{k}^{\text{opt}}}\end{split}$
(16)

where (b) is a direct application of Lemma 1. We define
$\mathbf{R}^{\text{sub}}=\left[R^{\text{sub}}_{1},\dots,R^{\text{sub}}_{K}\right]^{\intercal}$
and
$\mathbf{R}^{\text{opt}}=\left[R^{\text{opt}}_{1},\dots,R^{\text{opt}}_{K}\right]^{\intercal}$.
Due to the equivalence between (P1) and (P4), $\mathbf{R}^{\text{sub}}$ and
$\mathbf{R}^{\text{opt}}$ are, respectively, the approximate and optimal
solution to (P4). The series of inequalities in (15) and (16) mean that if
$\gamma$ is the approximation accuracy of (P2) in the sense of (12), then
$\gamma$ is also an approximation accuracy in the same sense for the
$\operatorname{SINR}$ and the rate at all locations. On the other hand, the
objective function of (P4) is the sum of strictly concave functions
$\log\left(e^{R_{k}}-1\right)$ with respect to $R_{k}$. These concave
functions are always upper bounded by the first order of their Taylor
expansion. Particularly, we have:

$\log\left(e^{R^{\text{sub}}_{k}}-1\right)\leq\log\left(e^{\gamma\lambda_{k}}-1\right)+\frac{e^{\gamma\lambda_{k}}}{e^{\gamma\lambda_{k}}-1}\left(R^{\text{sub}}_{k}-\gamma\lambda_{k}\right)$
(17)

for any arrival rate $\gamma\lambda_{k}$. Furthermore, for $x>0$, the function
$x\mapsto\log\left(e^{x}-1\right)$ is strictly increasing. This implies that
for any arrival rate vector $\bm{\lambda}$ inside the stability region, we can
write:

$\displaystyle\sum_{k=1}^{K}Y_{k}\log\left(e^{R^{\text{sub}}_{k}}-1\right)$
$\displaystyle\geq\sum_{k=1}^{K}Y_{k}\log\left(e^{\gamma\lambda_{k}}-1\right)$
(18)

because at all locations, $R^{\text{sub}}_{k}>\gamma\lambda_{k}$ is true for
all arrival rate vectors including $\mathbf{R}^{\text{opt}}$. Combining (17)
and (18), we obtain the following inequality:

$\sum_{k=1}^{K}Y_{k}\frac{e^{\gamma\lambda_{k}}}{e^{\gamma\lambda_{k}}-1}\left(\gamma\lambda_{k}-R^{\text{sub}}_{k}\right)\leq
0.$ (19)

We now introduce the quadratic Lyapunov function
$\mathcal{L}\left(\mathbf{Y}(t)\right)=\sum_{k=1}^{K}\frac{1}{2}\frac{e^{\gamma\lambda_{k}+\epsilon}}{e^{\gamma\lambda_{k}+\epsilon}-1}Y^{2}_{k}(t)$,
where for any arrival rate $\gamma\lambda_{k}$ lying strictly inside the
stability region, we select $\epsilon>0$ such that $\forall k$,
$\gamma\lambda_{k}+\epsilon$ is inside or on the boundary of the region. We
get from (14):

$\derivative{t}\mathcal{L}\left(\mathbf{Y}\right)=\sum_{k=1}^{K}Y_{k}\frac{e^{\gamma\lambda_{k}+\epsilon}}{e^{\gamma\lambda_{k}+\epsilon}-1}\left(\gamma\lambda_{k}-R^{\text{sub}}_{k}\right).$
(20)

Replacing $\gamma\lambda_{k}$ with $\gamma\lambda_{k}+\epsilon$ in (19)
yields:
$\sum_{k=1}^{K}Y_{k}\frac{e^{\gamma\lambda_{k}+\epsilon}}{e^{\gamma\lambda_{k}+\epsilon}-1}\left(\gamma\lambda_{k}+\epsilon-R^{\text{sub}}_{k}\right)\leq
0$. In view of (20), this reads:

$\derivative{t}\mathcal{L}\left(\mathbf{Y}\right)\leq-\epsilon\sum_{k=1}^{K}Y_{k}\frac{e^{\gamma\lambda_{k}+\epsilon}}{e^{\gamma\lambda_{k}+\epsilon}-1}.$
(21)

By using the the inequality

$2\times\mathcal{L}\left(\mathbf{Y}\right)=\sum_{k=1}^{K}\frac{e^{\gamma\lambda_{k}+\epsilon}}{e^{\gamma\lambda_{k}+\epsilon}-1}Y^{2}_{k}(t)\leq\left(\sum_{k=1}^{K}\frac{e^{\gamma\lambda_{k}+\epsilon}}{e^{\gamma\lambda_{k}+\epsilon}-1}Y_{k}(t)\right)^{2},$

we conclude that there exists a constant $\xi>0$ such that:

$\derivative{t}\mathcal{L}\left(\mathbf{Y}\right)\leq-\xi\sqrt{\mathcal{L}\left(\mathbf{Y}\right)}$
(22)

This implies for all $t\geq 0$ such that
$\mathcal{L}\left(\mathbf{Y}(t)\right)\geq 0$:

$\mathcal{L}\left(\mathbf{Y}(t)\right)\leq\left(\sqrt{\mathcal{L}\left(\mathbf{Y}(0)\right)}-\frac{\xi}{2}\,t\right)^{2}$
(23)

and also
$\sqrt{\mathcal{L}\left(\mathbf{Y}(t)\right)}\leq\left(\sqrt{\mathcal{L}\left(\mathbf{Y}(0)\right)}-\frac{\xi}{2}\,t\right)$
since $\mathcal{L}\left(\mathbf{Y}(t)\right)\geq 0$. Therefore, the inequality
in (22) means that if there exists some $T>0$ with
$\mathcal{L}\left(\mathbf{Y}(T)\right)=0$, then
$\mathcal{L}\left(\mathbf{Y}(t)\right)=0$ for all $t\geq T$. Looking at (23),
we select:

$T=\frac{2}{\xi}\sqrt{\sum_{k=1}^{K}\frac{1}{2}\frac{e^{\gamma\lambda_{k}+\epsilon}}{e^{\gamma\lambda_{k}+\epsilon}-1}}$
(24)

This implies that for all $t\geq T$, $\mathcal{L}\left(\mathbf{Y}(t)\right)$
and thus $\mathbf{Y}(t)$ are identically zero, and the system is stable. ∎

Figure 3: Simulation Results: Stability Region Figure 4: Simulation Results:
Moving average of the number of users for different arrival rates

Before providing our simulation results, we recapitulate the systematic way of
obtaining this result. First, we represent the queueing model by a vector
$\mathbf{X}(t)$ containing the number of users at each location. Given the
arrival model, we obtain the temporal equation describing the evolution of
$\mathbf{X}$. Then, using the law of large numbers, a similar limit equation
can be obtained for its fluid limit $\mathbf{Y}$. Knowing that the allocated
rate at each location should be strictly larger than the arrival rate in order
to guarantee the stability, we bound the suboptimal rates obtained from (P1)
and show that they are optimal in the scaled stability region. Finally, we use
the Lyapunov technique to show that $\mathbf{Y}$ reaches zero in finite time
for these arrival rates. It is then a direct conclusion that the network is
stable.

Figure 5: Simulation Results: Evolution of the stability metric for different
phase shift designs. The arrival rate is fixed at 0.1 users/slot at all
locations.

##  V NUMERICAL RESULTS

In this section, we report some simulation results in order to illustrate the
stability of the considered wireless system. We consider a $2\times 2$ km2
area with $N=128$ APs uniformly distributed in the lower left sub-region with
coordinates $x,y\in[-0.75,-1]$ and the RIS with $M=1600$ elements is located
at the origin. The bandwidth of the system is $20$ MHz and the carrier
frequency is $1.9$ GHz. The large scale fading coefficients are generated
according to the three-slope model in [27, Section VI].

For visual considerations, we start by considering $K=2$ locations in the
upper right region, with both coordinates being 0.25 and 0.75 for each
location respectively. The users arrive to each location according to a
Poisson process and transmit packets of $1$ Mbits with a power budget of $20$
dBm. In Figure 3 we compare the network stability region obtained by employing
our proposed phase shifts design and the one obtained by using time division
multiple access (TDMA) where at each timeslot, the RIS phases are optimized to
maximize the $\operatorname{SNR}$ of the active user. We notice that the TDMA
stability region is contained within the one guaranteed by our proposed RIS
configuration. In Figure 4, we show the moving average of the locations’ user
density. These correspond to two arrival vectors shown as points
$\bm{\lambda}_{1}$ and $\bm{\lambda}_{2}$ on the stability region in Figure 3.
Clearly, under the arrival rates of $\bm{\lambda}_{1}$ the network is stable
since the average traffic converges to a constant; while it is unstable for
the arrival rates of $\bm{\lambda}_{2}$ since this average diverges with time.
Thus $\bm{\lambda}_{2}$ is outside the stability region.

We now consider $K=400$ locations whose coordinates $x,y\in[0.05,1]$ form a
square mesh in the upper right sub-region. In Figure 5, we compare, for an
arrival rate equal to $0.1$ users/slot at all locations, the sum over the
locations of the moving average of the users density, when using our proposed
RIS configuration to the one obtained by drawing each phase uniformly at
random from $[0,2\pi]$ at each timeslot. It is clear that the network is not
stable using random RIS phases since the number of flows grows rapidly with
time, while this number is bounded when using the proposed phase shifts. In
Figure 6, we plot the stability metric computed at $T=10000$ timeslots versus
the arrival rate at all locations, for our proposed RIS phase shifts and the
random phases scheme. We first remark the metric varies exponentially with
respect to the arrival rate. Secondly, we evince that the network is stable
for a wider region of arrival rates when the RIS phases are optimized as we
proposed, compared to the case where they are randomly selected, since the
stability metric diverges at a prior point in the latter case.

Figure 6: Simulation Results: Stability metric versus the arrival rate at all
locations for different phase shift designs. The metric is calculated at
T=10000 timeslots.

##  VI CONCLUSION

In this paper, we studied the stability of a RIS-assisted cell-free MIMO
system where the number of users in the network is time-varying. We started by
obtaining a closed form physical layer rate expression, and then described the
data layer flow arrivals. Next, we proposed an optimization framework for the
RIS induced phases. With no power control nor scheduling, we showed that by
only optimizing the RIS phases, our proposed configuration ensures the
system’s stability inside a scaled version of the stability region. We have
also shown that this scaling factor is bounded.

## References

  * [1] T. L. Marzetta, “Noncooperative cellular wireless with unlimited numbers of base station antennas,” _IEEE transactions on wireless communications_ , vol. 9, no. 11, pp. 3590–3600, 2010.
  * [2] S. E. Hajri, M. Assaad, and G. Caire, “Scheduling in massive mimo: User clustering and pilot assignment,” in _2016 54th Annual Allerton Conference on Communication, Control, and Computing (Allerton)_. IEEE, 2016, pp. 107–114.
  * [3] F. Boccardi, R. W. Heath, A. Lozano, T. L. Marzetta, and P. Popovski, “Five disruptive technology directions for 5g,” _IEEE communications magazine_ , vol. 52, no. 2, pp. 74–80, 2014.
  * [4] C. Liaskos, S. Nie, A. Tsioliaridou, A. Pitsillides, S. Ioannidis, and I. Akyildiz, “A new wireless communication paradigm through software-controlled metasurfaces,” _IEEE Communications Magazine_ , vol. 56, no. 9, pp. 162–169, 2018.
  * [5] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 11, pp. 5394–5409, 2019.
  * [6] V. Croisfelt, F. Saggese, I. Leyva-Mayorga, R. Kotaba, G. Gradoni, and P. Popovski, “A random access protocol for ris-aided wireless communications,” in _2022 IEEE 23rd International Workshop on Signal Processing Advances in Wireless Communication (SPAWC)_. IEEE, 2022, pp. 1–5.
  * [7] X. Cao, B. Yang, C. Huang, G. C. Alexandropoulos, C. Yuen, Z. Han, H. V. Poor, and L. Hanzo, “Massive access of static and mobile users via reconfigurable intelligent surfaces: protocol design and performance analysis,” _IEEE Journal on Selected Areas in Communications_ , vol. 40, no. 4, pp. 1253–1269, 2022\.
  * [8] G. Yang, X. Xu, Y.-C. Liang, and M. Di Renzo, “Reconfigurable intelligent surface-assisted non-orthogonal multiple access,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 5, pp. 3137–3151, 2021.
  * [9] Z. Yang, J. Shi, Z. Li, M. Chen, W. Xu, and M. Shikh-Bahaei, “Energy efficient rate splitting multiple access (rsma) with reconfigurable intelligent surface,” in _2020 IEEE International Conference on Communications Workshops (ICC Workshops)_. IEEE, 2020, pp. 1–6.
  * [10] T. Van Chien, H. Q. Ngo, S. Chatzinotas, M. Di Renzo, and B. Ottersten, “Reconfigurable intelligent surface-assisted cell-free massive mimo systems over spatially-correlated channels,” _IEEE Transactions on Wireless Communications_ , vol. 21, no. 7, pp. 5106–5128, 2022.
  * [11] J. Denis and M. Assaad, “Improving cell-free massive mimo networks performance: A user scheduling approach,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 11, pp. 7360–7374, 2021.
  * [12] B. Al-Nahhas, M. Obeed, A. Chaaban, and M. J. Hossain, “Ris-aided cell-free massive mimo: Performance analysis and competitiveness,” in _IEEE International Conference on Communications Workshops (ICC Workshops)_. IEEE, 2021, pp. 1–6.
  * [13] M. Assaad, S. E. Hajri, T. Bonald, and A. Ephremides, “Power control in massive mimo with dynamic user population,” in _2018 IEEE Globecom Workshops (GC Wkshps)_. IEEE, 2018, pp. 1–6.
  * [14] E. Björnson and L. Sanguinetti, “Rayleigh fading modeling and channel hardening for reconfigurable intelligent surfaces,” _IEEE Wireless Communications Letters_ , vol. 10, no. 4, pp. 830–834, 2020.
  * [15] T. Van Chien, H. Q. Ngo, S. Chatzinotas, and B. Ottersten, “Reconfigurable intelligent surface-assisted massive mimo: Favorable propagation, channel hardening, and rank deficiency [lecture notes],” _IEEE Signal Processing Magazine_ , vol. 39, no. 3, pp. 97–104, 2022.
  * [16] Z. Wang, L. Liu, S. Zhang, and S. Cui, “Massive mimo communication with intelligent reflecting surface,” _arXiv preprint arXiv:2107.04255_ , 2021\.
  * [17] E. Björnson, J. Hoydis, L. Sanguinetti _et al._ , “Massive mimo networks: Spectral, energy, and hardware efficiency,” _Foundations and Trends® in Signal Processing_ , vol. 11, no. 3-4, pp. 154–655, 2017.
  * [18] M. J. Neely, “Stability and capacity regions or discrete time queueing networks,” _arXiv preprint arXiv:1003.3396_ , 2010.
  * [19] X. Zhang, _Matrix analysis and applications_. Cambridge University Press, 2017.
  * [20] J. Schur, “Bemerkungen zur theorie der beschränkten bilinearformen mit unendlich vielen veränderlichen.” _Journal für die reine und angewandte Mathematik_ , vol. 140, pp. 1–28, 1911. [Online]. Available: http://eudml.org/doc/149352
  * [21] A. M.-C. So, J. Zhang, and Y. Ye, “On approximating complex quadratic optimization problems via semidefinite programming relaxations,” _Mathematical Programming_ , vol. 110, no. 1, pp. 93–110, 2007.
  * [22] M. Grant and S. Boyd, “Cvx: Matlab software for disciplined convex programming, version 2.1,” 2014.
  * [23] Z.-Q. Luo, W.-K. Ma, A. M.-C. So, Y. Ye, and S. Zhang, “Semidefinite relaxation of quadratic optimization problems,” _IEEE Signal Processing Magazine_ , vol. 27, no. 3, pp. 20–34, 2010.
  * [24] T. Bonald and L. Massoulié, “Impact of fairness on internet performance,” in _Proceedings of the 2001 ACM SIGMETRICS international conference on Measurement and modeling of computer systems_ , 2001, pp. 82–91.
  * [25] M. Andrews, K. Kumaran, K. Ramanan, A. Stolyar, R. Vijayakumar, and P. Whiting, “Scheduling in a queuing system with asynchronously varying service rates,” _Probability in the Engineering and Informational Sciences_ , vol. 18, no. 2, pp. 191–217, 2004.
  * [26] P. S. Bullen, _Handbook of means and their inequalities_. Springer Science & Business Media, 2013, vol. 560.
  * [27] H. Q. Ngo, A. Ashikhmin, H. Yang, E. G. Larsson, and T. L. Marzetta, “Cell-free massive mimo versus small cells,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 3, pp. 1834–1850, 2017.