#  Robust Beamforming Design for Fairness-Aware Energy Efficiency Maximization
in RIS-Assisted mmWave Communications

Ahmed Magbool  Vaibhav Kumar
and Mark F. Flanagan  This publication has emanated from research conducted
with the financial support of Science Foundation Ireland under Grant Number
13/RC/2077_P2. Ahmed Magbool and Mark F. Flanagan are with the School of
Electrical and Electronic Engineering, University College Dublin, Belfield,
Dublin 4, Dublin 4, D04 V1W8, Ireland (e-mail<EMAIL_ADDRESS>mark.flanagan@ieee.org). Vaibhav Kumar was with the School of Electrical and
Electronic Engineering, University College Dublin, Dublin 4, D04 V1W8,
Ireland. He is now with Engineering Division, New York University Abu Dhabi
(NYUAD), Abu Dhabi 129188, UAE (e-mail: vaibhav.kumar@ieee.org).

###### Abstract

Users in millimeter-wave (mmWave) systems often exhibit diverse channel
strengths, which can negatively impact user fairness in resource allocation.
Moreover, exact channel state information (CSI) may not be available at the
transmitter, rendering suboptimal resource allocation. In this paper, we
address these issues within the context of energy efficiency maximization in
RIS-assisted mmWave systems. We first derive a tractable lower bound on the
achievable sum rate, taking into account CSI errors. Subsequently, we
formulate the optimization problem, targeting maximizing the system energy
efficiency while maintaining a minimum Jain’s fairness index controlled by a
tunable design parameter. The optimization problem is very challenging due to
the coupling of the optimization variables in the objective function and the
fairness constraint, as well as the existence of non-convex equality and
fractional constraints. To solve the optimization problem, we employ the
penalty dual decomposition method, together with a projected gradient ascent
based alternating optimization procedure. Simulation results demonstrate that
the proposed algorithm can achieve an optimal energy efficiency for a
prescribed Jain’s fairness index. In addition, adjusting the fairness design
parameter can yield a favorable trade-off between energy efficiency and user
fairness compared to methods that exclusively focus on optimizing one of these
metrics.

###### Index Terms:

Reconfigurable intelligent surfaces, mmWave communications, energy efficiency,
user fairness, imperfect CSI, projected gradient ascent, penalty dual
decomposition method.

##  I Introduction

The millimeter-wave (mmWave) band, extending from 28 GHz to 100 GHz, has been
exploited in existing fifth-generation (5G) wireless networks and is planned
for use in the upcoming sixth-generation (6G) wireless networks [1, 2].
However, signals at such high frequencies experience extreme propagation
conditions, limiting the system’s capabilities [3].

The utilization of reconfigurable intelligent surfaces (RISs) is regarded as a
promising solution to address the challenge of obscured line-of-sight (LoS) in
mmWave communications. RISs can establish a programmable propagation
environment by manipulating the phase of the incident signal through nearly
passive reflecting elements. Consequently, higher spectral efficiency can be
attained with reduced power consumption compared to the corresponding non-RIS-
aided systems [4, 5, 6].

Sustainability is a crucial aspect of 6G systems, aiming to achieve an energy
efficiency (EE) target of 1 Tbit/Joule [7, 8, 9]. With that in mind, the
escalating power consumption associated with massive data transmission within
6G networks has led researchers to focus on EE maximization for 6G networks
[10, 11]. A milestone work in EE optimization for RIS-aided systems was
presented in [12], offering two solutions to determine the optimal beamforming
design. The first solution used a gradient-based approach, while the second
relied on a fractional-programming-based alternating optimization (AO) method.
In [13], a genetic approach was proposed, employing a covariance matrix
adaptation evolution strategy and Dinkelbach’s method to maximize the system
EE for terahertz (THz) systems. Optimal beamformer design, coupled with an
asymptotic channel capacity characterization under hardware impairments, was
addressed in [14]. The problem of secrecy EE maximization in a multiple-input
multiple-output multiple eavesdropper (MIMOME) channel was tackled in [15],
where a stationary solution was obtained using a penalty-based approach. In
addition, a long short-term memory (LSTM)-based solution was explored to
design an energy-efficient RIS-aided system in [16].

However, the aforementioned works overlooked two crucial issues. The first is
the need to ensure an adequate degree of fairness among users. Since different
users may have distinct quality-of-service (QoS) requirements, proportional
fairness should be imposed, ensuring that the weighted rates of all users
together satisfy a fairness measure [17]. Ignoring this aspect leads to
allocating minimal resources to users with weaker channels, as their
contribution to the EE objective function is very minor [18]. This problem is
more pronounced in the mmWave band due to the severe path loss, resulting in
users close to the transmitter having much stronger channel gains than those
farther away from it. While incorporating rate constraints into the
optimization problem (e.g., by setting a specific threshold for each user’s
rate) can impose some level of fairness, simulation results in [12] indicate
that these threshold constraints can severely affect the EE. Moreover,
depending on the choice of these rate thresholds, a feasible solution to the
optimization problem may not exist. Second, while the above-referenced papers
assumed perfect channel state information (CSI) knowledge at the transmitter,
this assumption does not hold in practice.

In this paper, we target maximizing the system EE while ensuring user fairness
for mmWave multiuser systems with an array-of-subarrays (AoSA) hybrid
beamforming architecture. We choose this architecture because it is more
energy-efficient compared to the a fully-connected digital beamforming
architecture [19, 20]. The proposed method is robust to CSI imperfections and
relies only on the assumption of bounded CSI error. The main contributions of
this paper are summarized as follows:

  * •

For the first time in the literature, we present a solution to the problem of
EE maximization for a transmitter with AoSA architecture employing hybrid
beamforming in the context of RIS-assisted mmWave communications. Our solution
is sufficiently general to cover for the case of multiple-antenna receivers
capable of decoding multiple streams simultaneously. Unlike previous works on
this topic, which often assume a special precoder (e.g., the zero-forcing (ZF)
precoder in [12]), our approach addresses a more general case with the
presence of multi-user interference (MUI).

  * •

We address the practical scenario where perfect CSI is not accessible at the
transmitter. To model the channel estimation error, we adopt the widely-used
bounded-error model [21, 22], which constrains the estimation error of the
cascaded channel within a sphere of a specific radius. Employing this bounded
error model, we derive a tractable lower bound on the achievable sum rate
through the triangle and Cauchy–Schwarz inequalities.

  * •

We formulate the EE maximization problem with a constraint on users’ fairness
among users, characterized by Jain’s fairness index. The fairness constraint
imposes specific proportions between the users’ weighted rates, leading to an
improved system EE compared to the standard formulation found in the
literature, which typically involves minimum rate thresholds.

  * •

To address the challenging non-convex optimization problem, we employ a
penalty dual decomposition (PDD) method by incorporating the fairness
constraint into the objective function. Subsequently, we utilize a projected
gradient based AO algorithm to determine the optimal digital precoder, analog
precoder, RIS reflection matrix, and receive combiners that maximize the
system EE.

  * •

We present comprehensive numerical results to demonstrate the performance of
the proposed method. Our findings indicate that the proposed method can
achieve superior trade-offs between EE and fairness compared to the state-of-
the-art EE and fairness maximization techniques discussed in [12] and [23].

A related problem was previously addressed in our work [18]. However, that
study was limited to single-antenna receivers and MUI-free reception, while
the current work extends the scope to include the AoSA architecture, multiple-
antenna receivers, and the presence of MUI. Additionally, in [18], perfect CSI
knowledge is assumed to be available at the transmitter, whereas this work
emphasizes robust system design in the presence of channel estimation errors.
While [18] aimed at maximizing both EE and user fairness, this paper
concentrates on EE maximization with a fairness constraint. Given the nature
of the more complicated optimization problem in this work, we utilize the
projected gradient ascent and the PDD methods, which are fast and numerically
efficient approaches for large-scale systems. The proposed method does not
require specific convex solvers such as CVX or CVXPY.

The remainder of the paper is organized as follows. Section II presents the
system model. Section III formulates the EE maximization problem with fairness
constraint. Section IV presents a detailed description of the proposed
penalty-based AO framework. Simulation results and associated discussion are
presented in Section V. Finally, the paper is concluded in Section VI.

Notations: Bold lowercase and uppercase letters denote vectors and matrices,
respectively. $|\cdot|$ represents the magnitude of a complex number (for a
complex vector, it is assumed to operate element-wise). $[\mathbf{a}]_{i}$
denotes the $i$-th element of the vector $\mathbf{a}$.
$||\mathbf{\cdot}||_{2}$ and $||\mathbf{\cdot}||_{\mathsf{F}}$ represent the
Euclidean vector norm and the Frobenius matrix norm, respectively.
$\mathbf{(\cdot)}^{\mathsf{T}}$, $\mathbf{(\cdot)}^{\mathsf{H}}$ and
$\mathbf{(\cdot)}^{-1}$ denote the matrix transpose, matrix conjugate
transpose, and matrix inverse, respectively. $\mathbf{I}$ represents the
identity matrix. $\boldsymbol{1}$ denotes a column vector all of whose
elements are equal to one. $\text{diag}(a_{1},\dots,a_{N})$ denotes a diagonal
matrix with the elements $a_{1},\dots,a_{N}$ on the main diagonal and zeros
elsewhere, while $\text{blkdiag}(\mathbf{a}_{1},\dots,\mathbf{a}_{N})$
represents a block diagonal matrix with the vectors
$\mathbf{a}_{1},\dots,\mathbf{a}_{N}$ on the main block diagonal and zeros
elsewhere. $\mathbb{C}$ indicates the set of complex numbers, and
$j=\sqrt{-1}$ is the imaginary unit. $\mathbb{E}\\{\cdot\\}$ stands for the
expectation operator. $\mathcal{CN}(\mathbf{0},\mathbf{B})$ represents a
complex Gaussian random vector with a mean of $\mathbf{0}$ and a covariance
matrix of $\mathbf{B}$. $\nabla_{\mathbf{a}}f(\cdot)$ represents the complex-
valued gradient of the real-valued function $f(\cdot)$ with respect to
$\mathbf{a}$. $\odot$ and $\otimes$ denote the Khatri-Rao product and the
Kronecker product, respectively. $\mathcal{O}(\cdot)$ represents the big-O
notation to denote the computational complexity of an algorithm.

##  II System Model

We consider a multi-user downlink mmWave system in which a base-station (BS)
serves $K$ users as shown in Fig. 1. The BS needs to transmit $L$ data symbols
to each user, resulting in a total of $KL$ data streams. The BS employs hybrid
digital and analog precoding, where the analog precoder comprises $M$ radio-
frequency (RF) chains, each connected to a sub-array (SA) of $N_{\mathsf{T}}$
antenna elements (AEs). We assume that the direct paths between the BS and the
user equipment (UEs) are blocked and the communication is achieved via an RIS
with $N_{\text{RIS}}$ reflecting elements. At the receiver side, each UE is
assumed to have $N_{\mathsf{R}}$ AEs, and employs a digital combiner to decode
$L$ data streams, simultaneously111Due to constraints related to signal
processing capability, power consumption, and hardware design, the number of
streams for each UE (i.e., $L$) is typically small. In particular, the
following conditions must hold: $LK\leq M\leq N_{\mathsf{T}}$ and $L\leq
N_{\mathsf{R}}$.222This is a general system model and can be simplified into
several special cases. For instance, the case of single-antenna receivers can
be derived by setting $N_{\mathsf{R}}=L=1$, while $N_{\mathsf{T}}=1$
represents the scenario of fully digital antennas at the BS.. In the following
subsections, we describe the channel model, the signal model, and the channel
estimation error model.

Figure 1: System model for the proposed RIS-assisted mmWave system.

###  II-A Channel Model

We adopt the widely used far-field Saleh-Valenzuela (SV) channel model [24,
25, 26] to represent the BS-RIS and RIS-UE channels. We can write the BS-RIS
channel as

$\mathbf{H}_{\mathsf{T}}=\sum_{\ell=0}^{L_{\mathsf{T}}-1}\alpha_{\mathsf{T},\ell}\mathbf{a}_{\text{RIS}}(\vartheta_{\mathsf{R},\ell},\varphi_{\mathsf{R},\ell})\mathbf{a}_{\text{BS}}(\phi_{\mathsf{T},\ell})^{\mathsf{H}},$
(1)

where $L_{\mathsf{T}}$ is the number of paths between the BS and the RIS with
$\ell=0$ representing the LoS path and $\ell=1,\dots,L_{\mathsf{T}}-1$
representing the non-line-of-sight (NLoS) paths, and
$\alpha_{\mathsf{T},\ell}$ is the complex path gain of the $\ell$-th path. In
addition, $\mathbf{a}_{\text{BS}}(\phi_{\mathsf{T},\ell})$ is the beam
steering vector at the BS, which is a function of the angle of departure (AoD)
from BS $\phi_{\mathsf{T},\ell}$ of the $\ell$-th path from the BS. Assuming
that the distance between two adjacent AEs is $\lambda_{\mathrm{c}}/2$, with
$\lambda_{\mathrm{c}}$ indicating the carrier wavelength, and the distance
between the first antennas in two adjacent SAs is
$d_{\text{SA}}\lambda_{\mathrm{c}}/2$, the array response of the $n$-th AE at
the $m$-th SA at the BS can be written as

$a_{\text{BS}}^{m,n}(\phi_{\mathsf{T},\ell})=\exp\Big{(}{j\pi\big{(}(n-1)+(m-1)d_{\text{SA}}\big{)}\cos(\phi_{\mathsf{T},\ell}})\Big{)}.$
(2)

Also,
$\mathbf{a}_{\text{RIS}}(\vartheta_{\mathsf{R},\ell},\varphi_{\mathsf{R},\ell})$
is the beam steering vector of the RIS, which is a function of the azimuth and
elevation angles of arrival (AoA) to the RIS, $\vartheta_{\mathsf{R},\ell}$
and $\varphi_{\mathsf{R},\ell}$, respectively. Assuming that the horizontal
and vertical distances between two RIS elements is $\lambda_{\mathrm{c}}/2$,
the element in the $n_{x}$-th row and $n_{y}$-th column in the uniform planar
array (UPA) RIS have an array response of

$\begin{split}a_{\text{RIS}}^{n_{x},n_{y}}(\vartheta_{\mathsf{R},\ell},\varphi_{\mathsf{R},\ell})&=\exp\Big{(}j\pi\big{(}\varpi_{n_{x}}(\vartheta_{\mathsf{R},\ell},\varphi_{\mathsf{R},\ell})\\\
&+\varpi_{n_{y}}(\vartheta_{\mathsf{R},\ell},\varphi_{\mathsf{R},\ell})\big{)}\Big{)},\end{split}$
(3)

where
$\varpi_{n_{x}}(\vartheta_{\mathsf{R},\ell},\varphi_{\mathsf{R},\ell})=(n_{x}-1)\cos(\vartheta_{\mathsf{R},\ell})\sin(\varphi_{\mathsf{R},\ell})$
and
$\varpi_{n_{y}}(\vartheta_{\mathsf{R},\ell},\varphi_{\mathsf{R},\ell})=(n_{y}-1)\sin(\vartheta_{\mathsf{R},\ell})\sin(\varphi_{\mathsf{R},\ell})$.

Similarly, the channel between the RIS and the $k$-th UE is

$\mathbf{H}_{\mathsf{R},k}=\sum_{\ell=0}^{L_{\mathsf{R},k}-1}\alpha_{\mathsf{R},k,\ell}\mathbf{a}_{\text{UE},k}(\Xi_{\mathsf{R},k,\ell})\mathbf{a}_{\text{RIS}}(\vartheta_{\mathsf{T},\ell},\varphi_{\mathsf{T},\ell})^{\mathsf{H}},$
(4)

where $L_{\mathsf{R},k}$ is the number of paths between the RIS and the $k$-th
UE, $\alpha_{\mathsf{R},k,\ell}$ is the path gain of the $\ell$-th path
between the RIS and the $k$-th UE, and $\vartheta_{\mathsf{T},\ell}$ and
$\varphi_{\mathsf{T},\ell}$ are the azimuth and elevation AoD from the RIS,
respectively. Additionally,
$\mathbf{a}_{\text{UE},k}(\Xi_{\mathsf{R},k,\ell})$ is array response vector
for the $k$-th UE, which can be written as (assuming $\lambda_{\mathrm{c}}/2$
spacing between AEs)

$\mathbf{a}_{\text{UE},k}(\Xi_{\mathsf{R},k,\ell})=\big{[}1,\dots,\exp\big{(}j\pi(N_{\mathsf{R}}-1)\cos(\Xi_{\mathsf{R},k,\ell})\big{)}\big{]},$
(5)

where $\Xi_{\mathsf{R},k,\ell}$ is the AoA of the $\ell$-th path to the $k$-th
UE.

###  II-B Signal Model

The signal transmitted by the BS can be expressed as

$\mathbf{x}=\mathbf{A}\mathbf{D}\mathbf{s},$ (6)

where $\mathbf{s}=[s_{1,1},\dots,s_{K,L}]^{\mathsf{T}}\in\mathbb{C}^{KL\times
1}$ is a vector containing the transmitted symbols, with $s_{k,\ell}$
representing the $\ell$-th symbol to be transmitted to the $k$-th user, where
$k\in\\{1,2,...,K\\}$ and $\ell\in\\{1,2,...,L\\}$, and
$\mathbb{E}\\{\mathbf{s}\mathbf{s}^{\mathsf{H}}\\}=\mathbf{I}$. Also,
$\mathbf{D}=[\mathbf{d}_{1,1},\dots,\mathbf{d}_{K,L}]\in\mathbb{C}^{M\times
KL}$ is the digital precoder satisfying the power budget constraint, i.e.,
$||\mathbf{D}||_{\mathsf{F}}^{2}\leq P_{\text{max}}$, with $P_{\text{max}}$
denoting the total transmit power budget and
$\mathbf{d}_{k,\ell}\in\mathbb{C}^{M\times 1}$ denoting the digital precoder
for the $\ell$-th data stream intended for the $k$-th user. In addition,
$\mathbf{A}=\text{blkdiag}(\mathbf{a}_{1},\dots,\mathbf{a}_{M})\in\mathbb{C}^{MN_{\mathsf{T}}\times
M}$ is the analog precoder, with
$\mathbf{a}_{m}\in\mathbb{C}^{N_{\mathsf{T}}\times 1}$ representing the analog
precoder employed by the $m$-th RF chain, with elements obeying the constant
modulus (CM) constraint, i.e.,
$|[\mathbf{a}_{m}]_{i}|=1/\sqrt{N_{\mathsf{T}}}$.

At the other end, the post-combining received signal by the $k$-th user to
decode the $\ell$-th data stream is

$y_{k,\ell}=\mathbf{c}_{k,\ell}^{\mathsf{H}}\mathbf{G}_{k}(\boldsymbol{\Theta})\mathbf{x}+\mathbf{c}_{k,\ell}^{\mathsf{H}}\mathbf{n}_{k},$
(7)

where $\mathbf{c}_{k,\ell}\in\mathbb{C}^{N_{\mathsf{R}}\times 1}$ is the
combiner employed by the $k$-th user to decode the $\ell$-th data stream that
satisfies the receive power budget constraint
$\sum_{\ell=1}^{L}||\mathbf{c}_{k,\ell}||_{2}^{2}=P_{\mathsf{R},k}$, and
$\mathbf{n}_{k}\sim\mathcal{CN}(\mathbf{0},\sigma^{2}\mathbf{I})$ is the
complex additive white Gaussian noise (AWGN). Furthermore,
$\mathbf{G}_{k}(\boldsymbol{\Theta})\in\mathbb{C}^{N_{\mathsf{R}}\times
MN_{\mathsf{T}}}$ denotes the overall channel matrix for the $k$-th user,
which can be written as

$\mathbf{G}_{k}(\boldsymbol{\Theta})=\mathbf{H}_{\mathsf{R},k}\boldsymbol{\Theta}\mathbf{H}_{\mathsf{T}},$
(8)

where
$\boldsymbol{\Theta}=\text{diag}(\theta_{1},\dots,\theta_{N_{\text{RIS}}})\in\mathbb{C}^{N_{\text{RIS}}\times
N_{\text{RIS}}}$ is the RIS reflection coefficient with $|\theta_{n}|=1$ for
all $n\in\\{1,\dots,N_{\text{RIS}}\\}$. However, in practice, due to channel
estimation errors, perfect CSI will in general not be available at the
transmitter. In the next subsection, we introduce the channel estimation error
model.

###  II-C Channel Estimation Error Model

To characterize the channel estimation error model, we express the channel of
the $k$-th user as

$\mathbf{G}_{k}(\boldsymbol{\Theta})=\hat{\mathbf{G}}_{k}(\boldsymbol{\Theta})+\boldsymbol{\Lambda}_{k}(\boldsymbol{\Theta}),$
(9)

where $\hat{\mathbf{G}}_{k}(\boldsymbol{\Theta})$ represents the estimated
channel matrix at the BS, and $\boldsymbol{\Lambda}_{k}(\boldsymbol{\Theta})$
denotes the error matrix resulting from imperfect estimation of
$\mathbf{G}_{k}(\boldsymbol{\Theta})$.

As RISs are inherently passive, they lack the capability of directly
transmitting or processing pilots. Consequently, most channel estimation
algorithms focus on estimating the cascaded channel [27], which can be
obtained by vectorizing (9) as

$\text{vec}(\mathbf{G}_{k}(\boldsymbol{\Theta}))=\underbrace{(\mathbf{H}_{\mathsf{T}}^{\mathsf{T}}\odot\mathbf{H}_{\mathsf{R},k})}_{\mathbf{H}_{k}}\boldsymbol{\theta}=(\hat{\mathbf{H}}_{k}+\boldsymbol{\Delta}_{k})\boldsymbol{\theta},$
(10)

where $\boldsymbol{\theta}=\text{diag}(\boldsymbol{\Theta})$,
$\hat{\mathbf{H}}_{k}$ is the estimated version of $\mathbf{H}_{k}$, and
$\boldsymbol{\Delta}_{k}$ is the error matrix incurred during the estimation
of $\mathbf{H}_{k}$. The relationship between $\boldsymbol{\Delta}_{k}$ and
$\boldsymbol{\Lambda}_{k}(\boldsymbol{\Theta})$ in (9) can be expressed as
$\text{vec}\big{(}\boldsymbol{\Lambda}_{k}(\boldsymbol{\Theta})\big{)}=\boldsymbol{\Delta}_{k}\boldsymbol{\theta}$.

In this paper, we assume that the channel estimation error is bounded within a
sphere of radius $\delta_{k}$. This bounded error model is widely adopted in
the literature due to its effectiveness in characterizing quantization errors
which lie within a bounded region [21, 22].

##  III Problem Formulation

In this section, we formulate the problem of maximizing the system EE while
adhering to a fairness constraint. We begin by expressing the rate of the
$\ell$-th stream decoded by the $k$-th user as

$\begin{split}&R_{k,\ell}(\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\mathbf{c}_{k,\ell})=\log_{2}\bigg{(}1+\\\
&\frac{|\mathbf{c}_{k,\ell}^{\mathsf{H}}\mathbf{G}_{k}(\boldsymbol{\Theta})\mathbf{A}\mathbf{d}_{k,\ell}|^{2}}{\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{{v=1}}^{L}}|\mathbf{c}_{k,\ell}^{\mathsf{H}}\mathbf{G}_{k}(\boldsymbol{\Theta})\mathbf{A}\mathbf{d}_{i,v}|^{2}+\mathbb{E}\\{|\mathbf{c}_{k,\ell}^{\mathsf{H}}\mathbf{n}_{k}|^{2}\\}}\bigg{)}.\end{split}$
(11)

In addition, the total power consumption of the system can be expressed as
(c.f. [12])

$\begin{split}P_{\text{tot}}(\mathbf{D})&=\underbrace{P_{\text{BS}}+\xi||\mathbf{D}||_{\mathsf{F}}^{2}+MN_{\mathsf{T}}P_{\text{RF,T}}}_{\text{BS
power consumption}}\\\ &+\underbrace{N_{\text{RIS}}P_{\theta}}_{\text{RIS
power
consumption}}+\underbrace{\sum_{k=1}^{K}(P_{\text{UE},k}+P_{\text{R},k})}_{\text{UEs
power consumption}},\end{split}$ (12)

where $P_{\text{BS}}$ and $P_{\text{UE},k}$ are the total hardware static
power consumption at the BS and at the $k$-th UE, respectively, $\xi$ is the
power amplification factor at the BS, $P_{\text{RF,T}}$ and $P_{\theta}$ are
the power consumption of each phase shifter at the BS and at the RIS,
respectively.

We can then express the system EE as the sum rate divided by the total power
consumption, i.e.,

$\eta(\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\mathbf{C})=\frac{R_{\text{sum}}(\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\mathbf{C})}{P_{\text{tot}}(\mathbf{D})},$
(13)

where
$R_{\text{sum}}(\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\mathbf{C})\triangleq\sum_{K=1}^{K}\sum_{\ell=1}^{L}R_{k,\ell}(\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\mathbf{c}_{k,\ell})$
and $\mathbf{C}\triangleq[\mathbf{c}_{1,1},\dots,\mathbf{c}_{K,L}]$. The goal
is to allocate the available resources to maximize the system EE. However, as
the channel estimation error matrices
$\boldsymbol{\Delta}_{1},\dots,\boldsymbol{\Delta}_{K}$ are not accessible at
the transmitter, we introduce the following proposition.

######  Proposition 1.

The rate of the $\ell$-th data stream decoded by the $k$-th user can be lower
bounded as

$\begin{split}&R_{k,\ell}(\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\mathbf{c}_{k,\ell})\geq\underline{R}_{k,\ell}(\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\tilde{\mathbf{c}}_{k,\ell})\triangleq\log_{2}\bigg{(}1+\\\
&\frac{\varrho_{k}|\tilde{\mathbf{c}}_{k,\ell}^{\mathsf{H}}\hat{\mathbf{G}}_{k}(\boldsymbol{\Theta})\mathbf{A}\mathbf{d}_{k,\ell}|^{2}}{\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}\big{(}|\tilde{\mathbf{c}}_{k,\ell}^{\mathsf{H}}\hat{\mathbf{G}}_{k}(\boldsymbol{\Theta})\mathbf{A}\mathbf{d}_{i,v}|+\delta_{k}\sqrt{N_{\mathsf{RIS}}}||\mathbf{d}_{i,v}||_{2}\big{)}^{2}+\sigma^{2}}\bigg{)},\end{split}$
(14)

where
$\varrho_{k}\triangleq\big{(}1-\frac{\delta_{k}}{||\hat{\mathbf{H}}_{k}||_{\mathsf{F}}}\big{)}^{2}$
and
$\tilde{\mathbf{c}}_{k,\ell}\triangleq\mathbf{c}_{k,\ell}/||\mathbf{c}_{k,\ell}||_{2}$.

###### Proof.

See Appendix A. ∎

Using Proposition 1, we can obtain a lower bound on the system EE as

$\eta(\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\mathbf{C})\geq\underline{\eta}(\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\tilde{\mathbf{C}})\triangleq\frac{\underline{R}_{\text{sum}}(\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\tilde{\mathbf{C}})}{P_{\text{tot}}(\mathbf{D})},$
(15)

where
$\underline{R}_{\text{sum}}(\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\tilde{\mathbf{C}})\triangleq\sum_{k=1}^{K}\sum_{\ell=1}^{L}\underline{R}_{k,\ell}(\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\tilde{\mathbf{c}}_{k,\ell})$
and
$\tilde{\mathbf{C}}\triangleq[\tilde{\mathbf{c}}_{1,1},\dots,\tilde{\mathbf{c}}_{K,L}]$.

While it is crucial from a system design perspective to maintain high EE by
maximizing (15), ensuring acceptable QoS remains a primary concern for end-
users. However, users in mmWave communication systems often exhibit
considerable variations in their channel strengths due to the severe path loss
[3]. In such cases, weaker users, unable to significantly contribute to the EE
objective function (15), "are likely to be treated as unwanted interference to
stronger users, resulting in a limited resource allocation for the weaker
ones. To address this issue, we aim to maximize the system EE while ensuring a
certain level of fairness among users. This objective is pursued by tackling
the following optimization problem:

$\displaystyle\max_{\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\tilde{\mathbf{C}}}\
\ \
\underline{\eta}(\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\tilde{\mathbf{C}})$
(16a) $\displaystyle\text{s.t.}\ \
\rho\leq\frac{\big{(}\sum_{k=1}^{K}\frac{1}{w_{k}}\sum_{\ell=1}^{L}\underline{R}_{k,\ell}(\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\tilde{\mathbf{c}}_{k,\ell})\big{)}^{2}}{K\sum_{k=1}^{K}\big{(}\frac{1}{w_{k}}\sum_{\ell=1}^{L}\underline{R}_{k,\ell}(\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\tilde{\mathbf{c}}_{k,\ell})\big{)}^{2}}$
(16b) $\displaystyle||\mathbf{D}||_{\mathsf{F}}^{2}\leq P_{\text{max}},$ (16c)
$\displaystyle|\theta_{n}|=1,\ \forall n\in\\{1,\dots,N_{\text{RIS}}\\},$
(16d) $\displaystyle|[\mathbf{a}_{m}]_{i}|=\frac{1}{\sqrt{N_{\mathsf{T}}}},\
\forall m\in\\{1,\dots,M\\},\ \forall i\in\\{1,\dots,N_{\mathsf{T}}\\},$ (16e)
$\displaystyle||\tilde{\mathbf{c}}_{k,\ell}||_{2}=1,\ \forall
k\in\\{1,\dots,K\\},\ \forall\ell\in\\{1,\dots,L\\},$ (16f)

where $w_{k}$ is a weight assigned to the $k$-th user by the system provider,
representing the QoS for that user. Additionally, $\rho\in[1/K,1]$ serves as a
design parameter incorporated to manage the Jain’s fairness index of the
system. Jain’s fairness index is a widely recognized metric that quantifies
the fairness level among proportional rates in a system, ranging between $1/K$
and 1. A Jain’s fairness index of 1 indicates perfect fairness, implying that
all users have the same proportional rate [28].

Solving (16) poses a significant challenge due to the coupling between the
optimization variables in the objective function and the fairness constraint.
Moreover, handling the non-convex equality constraints (16d), (16e) and (16f)
is challenging. In the following section, we introduce projected gradient
ascent based AO procedures to address this optimization problem.

##  IV Proposed Solution

In this section, we present a penalty-based method to get rid of the
complicated fairness constraint, which otherwise is very challenging to handle
due to the fractional form, and also due to the coupling between the design
variables. We then apply an AO-based approach to update each of the design
variables in a one-by-one fashion.

The general framework for updating each variable comprises three main steps:

1) Obtaining the augmented Lagrangian function: First, to address the
challenging fairness constraint (16b), the PDD method [29] is employed. The
PDD method operates by augmenting challenging constraints with penalty into
the objective function, and ensures that both the original and the transformed
problems share the same stationary solution. For the problem at hand, to
obtain the penalty function, we first rewrite the constraint (16b) as

$\rho
K\sum_{k=1}^{K}\Big{(}\frac{1}{w_{k}}\sum_{\ell=1}^{L}\underline{R}_{k,\ell}(\boldsymbol{\chi})\Big{)}^{2}-\Big{(}\sum_{k=1}^{K}\frac{1}{w_{k}}\sum_{\ell=1}^{L}\underline{R}_{k,\ell}(\boldsymbol{\chi})\Big{)}^{2}\leq
0,$ (17)

where
$\boldsymbol{\chi}\in\\{\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\tilde{\mathbf{C}}\\}$
is the relevant optimization variable333In this paper, having some arguments
omitted in a multi-variable function indicates that their values are being
held fixed.. Based on (17), we define the following penalty function:

$\begin{split}\mathcal{G}(\boldsymbol{\chi},\mu)&\triangleq\rho
K\sum_{k=1}^{K}\Big{(}\frac{1}{w_{k}}\sum_{\ell=1}^{L}\underline{R}_{k,\ell}(\boldsymbol{\chi})\Big{)}^{2}\\\
&-\Big{(}\sum_{k=1}^{K}\frac{1}{w_{k}}\sum_{\ell=1}^{L}\underline{R}_{k,\ell}(\boldsymbol{\chi})\Big{)}^{2}+\mu,\end{split}$
(18)

where $\mu\geq 0$ is a slack variable used to eliminate the penalty function
whenever the condition (17) is satisfied.

Then, the PDD method converts the following constrained optimization problem:

$\displaystyle\max_{\boldsymbol{\chi}}\ \ \ $
$\displaystyle\underline{\eta}(\boldsymbol{\chi})$ (19a) s.t.
$\displaystyle\eqref{eq:con_re},$ (19b)

into the following form:

$\displaystyle\max_{\boldsymbol{\chi},\mu\geq 0}\ \ \ $
$\displaystyle\mathcal{H}_{\gamma,\omega}(\boldsymbol{\chi},\mu),$ (20a)

where $\mathcal{H}_{\gamma,\omega}(\boldsymbol{\chi},\mu)$ is the augmented
Lagrangian function in the form

$\mathcal{H}_{\gamma,\omega}(\boldsymbol{\chi},\mu)\triangleq\underline{\eta}(\boldsymbol{\chi})-\gamma\mathcal{G}(\boldsymbol{\chi},\mu)-\frac{1}{2\omega}\mathcal{G}^{2}(\boldsymbol{\chi},\mu),$
(21)

where $\gamma$ is the Lagrange multiplier and $\omega>0$ is a penalty
parameter [29].

For fixed values of $\gamma$ and $\omega$, the optimal solution of the
optimization problem (20) is obtained using a projection-based AO procedure,
as will be demonstrated later in this section. Once convergence is achieved
for specific values of $\gamma$ and $\omega$, they can be updated as [29]:

$\gamma\leftarrow\gamma+\frac{1}{\omega}\mathcal{H}_{\gamma,\omega}(\mathbf{D},\mathbf{A},\boldsymbol{\Theta},\tilde{\mathbf{C}},\mu),$
(22)

and

$\omega\leftarrow\psi\omega,$ (23)

where $\psi\in(0,1]$ is a design parameter that controls the speed of
attenuation of the constraint violation.

2) Gradient ascent update: To update the relevant variable
$\boldsymbol{\chi}$, the gradient of the augmented Lagrangian function with
respect to the optimization variable (i.e.,
$\nabla_{\boldsymbol{\chi}}\mathcal{H}_{\gamma,\omega}(\boldsymbol{\chi})$) is
obtained. Then, the variable is updated via

$\bar{\boldsymbol{\chi}}^{(t+1)}=\boldsymbol{\chi}^{(t)}+\alpha_{\boldsymbol{\chi}}^{(t)}\nabla_{\boldsymbol{\chi}}\mathcal{H}_{\gamma,\omega}(\boldsymbol{\chi}^{(t)}),$
(24)

where $t$ is the iteration number, $\alpha_{\boldsymbol{\chi}}^{(t)}>0$ is the
step size for updating $\boldsymbol{\chi}$ at the $t$-th iteration, and
$\bar{\boldsymbol{\chi}}^{(t+1)}$ and $\boldsymbol{\chi}^{(t)}$ are the values
of $\boldsymbol{\chi}$ at the $(t+1)$-th and the $t$-th iterations of the
gradient ascent algorithm, respectively.

3) Projection onto the feasible set: The updated value of the variable
$\boldsymbol{\chi}$, i.e., $\bar{\boldsymbol{\chi}}^{(t+1)}$ may not lie in
its feasible set. Therefore, we project the updated variable orthogonally onto
the feasible set using the projection function:

$\boldsymbol{\chi}^{(t+1)}=\mathcal{T}_{\boldsymbol{\chi}}(\bar{\boldsymbol{\chi}}^{(t+1)}).$
(25)

We detail the procedure of updating each variable in the following
subsections.

###  IV-A Digital Precoder

We first assume that all variables other than the digital precoder are fixed
(i.e., the values of $\mathbf{A},\boldsymbol{\Theta},\tilde{\mathbf{C}}$ and
$\mu$ are held fixed). We can then write the EE maximization problem with
respect to the digital precoder as

$\displaystyle\max_{\mathbf{D}}\ \ \ $
$\displaystyle\underline{\eta}(\mathbf{D})$ (26a) s.t.
$\displaystyle\eqref{F_const},\ \eqref{Pow_con}.\vspace{-0.2cm}$ (26b)

The augmented Lagrangian function after incorporating the constraint (16b)
into (26a) is

$\mathcal{H}_{\gamma,\omega}(\mathbf{D})=\underline{\eta}(\mathbf{D})-\gamma\mathcal{G}(\mathbf{D})-\frac{1}{2\omega}\mathcal{G}^{2}(\mathbf{D}),$
(27)

To find the gradient of (27) with repect to $\mathbf{D}$, we begin by re-
expressing $\underline{R}_{k,\ell}(\mathbf{D})$ in (14) as

$\begin{split}&\underline{R}_{k,\ell}(\mathbf{D})=\log_{2}\Big{(}\varrho_{k}|\mathbf{b}_{k,\ell}^{\mathsf{H}}\mathbf{d}_{k,\ell}|^{2}\\\
&+\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}\big{(}|\mathbf{b}_{k,\ell}^{\mathsf{H}}\mathbf{d}_{i,v}|+\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2}\big{)}^{2}+\sigma^{2}\Big{)}\\\
&-\log_{2}\Big{(}\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}\big{(}|\mathbf{b}_{k,\ell}^{\mathsf{H}}\mathbf{d}_{i,v}|+\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2}\big{)}^{2}+\sigma^{2}\Big{)},\end{split}$
(28)

where
$\mathbf{b}_{k,\ell}=\mathbf{A}^{\mathsf{H}}\hat{\mathbf{G}}_{k}^{\mathsf{H}}(\boldsymbol{\Theta})\tilde{\mathbf{c}}_{k,\ell}$.

We can obtain the gradient of $\underline{R}_{k,\ell}(\mathbf{d}_{m,u})$ with
respect to $\mathbf{d}_{m,u}$ by considering two cases. First, when
$[m,u]=[k,\ell]$, we have the gradient expression shown in (29).

$\nabla_{\mathbf{d}_{m,u}}\underline{R}_{k,\ell}(\mathbf{D})\Big{|}_{[m,u]=[k,\ell]}=\frac{2}{\ln
2}\
\frac{\varrho_{k}\mathbf{b}_{k,\ell}\mathbf{b}_{k,\ell}^{\mathsf{H}}\mathbf{d}_{k,\ell}}{\varrho_{k}|\mathbf{b}_{k,\ell}^{\mathsf{H}}\mathbf{d}_{k,\ell}|^{2}+\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}(|\mathbf{b}_{k,\ell}^{\mathsf{H}}\mathbf{d}_{i,v}|+\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2})^{2}+\sigma^{2}}.$
(29)

On the other hands, when $[m,u]\neq[k,\ell]$, we can obtain the gradient
expression shown in (30).

$\begin{split}\nabla_{\mathbf{d}_{m,u}}\underline{R}_{k,\ell}(\mathbf{D})\Big{|}_{[m,u]\neq[k,\ell]}=\frac{2}{\ln
2}\bigg{[}&\frac{\big{(}1+\frac{\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{m,u}||_{2}}{|\mathbf{b}_{k,\ell}^{\mathsf{H}}\mathbf{d}_{m,u}|}\big{)}\mathbf{b}_{k,\ell}\mathbf{b}_{k,\ell}^{\mathsf{H}}+\big{(}\delta_{k}^{2}N_{\text{RIS}}+\frac{\delta_{k}\sqrt{N_{\text{RIS}}}|\mathbf{b}_{k,\ell}^{\mathsf{H}}\mathbf{d}_{m,u}|}{||\mathbf{d}_{m,u}||_{2}}\big{)}\mathbf{I}}{\varrho_{k}|\mathbf{b}_{k,\ell}^{\mathsf{H}}\mathbf{d}_{k,\ell}|^{2}+\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}(|\mathbf{b}_{k,\ell}^{\mathsf{H}}\mathbf{d}_{i,v}|+\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2})^{2}+\sigma^{2}}\\\
&-\frac{\big{(}1+\frac{\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{m,u}||_{2}}{|\mathbf{b}_{k,\ell}^{\mathsf{H}}\mathbf{d}_{m,u}|}\big{)}\mathbf{b}_{k,\ell}\mathbf{b}_{k,\ell}^{\mathsf{H}}+\big{(}\delta_{k}^{2}N_{\text{RIS}}+\frac{\delta_{k}\sqrt{N_{\text{RIS}}}|\mathbf{b}_{k,\ell}^{\mathsf{H}}\mathbf{d}_{m,u}|}{||\mathbf{d}_{m,u}||_{2}}\big{)}\mathbf{I}}{\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}(|\mathbf{b}_{k,\ell}^{\mathsf{H}}\mathbf{d}_{i,v}|+\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2})^{2}+\sigma^{2}}\bigg{]}\mathbf{d}_{m,u}.\end{split}$
(30)

Thus, the gradient of $\underline{R}_{k,\ell}(\mathbf{D})$ with respect to
$\mathbf{d}_{m,u}$ can be expressed as

$\nabla_{\mathbf{d}_{m,u}}\underline{R}_{k,\ell}(\mathbf{D})=\begin{cases}\eqref{eq:grad_c1},&\text{if}\
[m,u]=[k,\ell],\\\ \eqref{eq:grad_c2},&\text{if}\
[m,u]\neq[k,\ell].\end{cases}$ (31)

The gradient of $\underline{R}_{k,\ell}(\mathbf{D})$ with respect to
$\mathbf{D}$ can then be written as

$\nabla_{\mathbf{D}}\underline{R}_{k,\ell}(\mathbf{D})=\big{[}\nabla_{\mathbf{d}_{1,1}}\underline{R}_{k,\ell}(\mathbf{D}),\dots,\nabla_{\mathbf{d}_{K,L}}\underline{R}_{k,\ell}(\mathbf{D})\big{]},$
(32)

and that of $\underline{R}_{\text{sum}}(\mathbf{D})$ with respect to
$\mathbf{D}$ is

$\nabla_{\mathbf{D}}\underline{R}_{\text{sum}}(\mathbf{D})=\sum_{k=1}^{K}\sum_{\ell=1}^{L}\nabla_{\mathbf{D}}\underline{R}_{k,\ell}(\mathbf{D}).$
(33)

We can also obtain the gradient of $P_{\text{tot}}(\mathbf{D})$ with respect
to $\mathbf{D}$ as

$\nabla_{\mathbf{D}}P_{\text{tot}}(\mathbf{D})=2\xi\mathbf{D}.$ (34)

Then the gradient of $\underline{\eta}(\mathbf{D})$ with respect to
$\mathbf{D}$ is given by

$\nabla_{\mathbf{D}}\underline{\eta}(\mathbf{D})=\frac{P_{\text{tot}}(\mathbf{D})\nabla_{\mathbf{D}}\underline{R}_{\text{sum}}(\mathbf{D})-\underline{R}_{\text{sum}}(\mathbf{D})\nabla_{\mathbf{D}}P_{\text{tot}}(\mathbf{D})}{P_{\text{tot}}^{2}(\mathbf{D})}.$
(35)

We also need to find the gradient of $\mathcal{G}(\mathbf{D})$ with respect to
$\mathbf{D}$, which can be written as

$\begin{split}&\nabla_{\mathbf{D}}\mathcal{G}(\mathbf{D})=2\rho
K\Big{(}\sum_{k=1}^{K}\frac{1}{w_{k}^{2}}\sum_{\ell=1}^{L}\underline{R}_{k,\ell}(\mathbf{D})\sum_{m=1}^{L}\nabla_{\mathbf{D}}\underline{R}_{k,m}(\mathbf{D})\Big{)}\\\
&-2\Big{(}\sum_{k=1}^{K}\frac{1}{w_{k}}\sum_{\ell=1}^{L}\underline{R}_{k,\ell}(\mathbf{D})\Big{)}\Big{(}\sum_{i=1}^{K}\frac{1}{w_{i}}\sum_{m=1}^{L}\nabla_{\mathbf{D}}\underline{R}_{i,m}(\mathbf{D})\Big{)}.\end{split}$
(36)

Thus, the gradient of $\mathcal{H}_{\gamma,\omega}(\mathbf{D})$ with respect
to $\mathbf{D}$ is

$\nabla_{\mathbf{D}}\mathcal{H}_{\gamma,\omega}(\mathbf{D})=\nabla_{\mathbf{D}}\underline{\eta}(\mathbf{D})-\gamma\nabla_{\mathbf{D}}\mathcal{G}(\mathbf{D})-\frac{1}{\omega}\mathcal{G}(\mathbf{D})\nabla_{\mathbf{D}}\mathcal{G}(\mathbf{D}),$
(37)

and the update of $\mathbf{D}$ should follow

$\bar{\mathbf{D}}^{(t+1)}=\mathbf{D}^{(t)}+\alpha_{\mathbf{D}}^{(t)}\nabla_{\mathbf{D}}\mathcal{H}_{\gamma,\omega}(\mathbf{D}^{(t)}).$
(38)

Next, we note that the projection of the updated variable
$\bar{\mathbf{D}}^{(t+1)}$ onto the constraint (16c) is given by

$\mathbf{D}^{(t+1)}=\min\bigg{(}1,\frac{\sqrt{P_{\text{max}}}}{||\bar{\mathbf{D}}^{(t+1)}||_{\mathsf{F}}}\bigg{)}\bar{\mathbf{D}}^{(t+1)}.$
(39)

###  IV-B Analog Precoder

Assuming that the values of $\mathbf{D},\boldsymbol{\Theta}$ and
$\tilde{\mathbf{C}}$ are fixed, the EE optimization for the analog precoder is

$\displaystyle\max_{\mathbf{A}}\ \ \ $
$\displaystyle\underline{\eta}(\mathbf{A}),$ (40a) s.t.
$\displaystyle\eqref{F_const},\ \eqref{CM_BS}.$ (40b)

We first express the augmented Lagrangian function of (40a) when augmenting
the constraint (16b) into the objective function as

$\mathcal{H}_{\gamma,\omega}(\mathbf{A})=\frac{1}{P_{\text{tot}}}\underline{R}_{\text{sum}}(\mathbf{A})-\gamma\mathcal{G}(\mathbf{A})-\frac{1}{2\omega}\mathcal{G}^{2}(\mathbf{A}).$
(41)

To derive the gradient expression of the sum rate with respect to
$\mathbf{A}$, we begin by expressing the $k$-th user’s hybrid precoding vector
for the $\ell$-th stream as

$\mathbf{f}_{k,\ell}=\mathbf{A}\mathbf{d}_{k,\ell}=\big{[}[\mathbf{d}_{k,\ell}]_{1}\mathbf{a}_{1}^{\mathsf{T}},\dots,[\mathbf{d}_{k,\ell}]_{M}\mathbf{a}_{M}^{\mathsf{T}}\big{]}^{\mathsf{T}}=\tilde{\mathbf{D}}_{k,\ell}\mathbf{a},$
(42)

where
$\tilde{\mathbf{D}}_{k,\ell}=\text{diag}(\big{[}[\mathbf{d}_{k,\ell}]_{1}\boldsymbol{1}^{\mathsf{T}},\dots,[\mathbf{d}_{k,\ell}]_{M}\boldsymbol{1}^{\mathsf{T}}\big{]}^{\mathsf{T}})$,
and
$\mathbf{a}=[\mathbf{a}_{1}^{\mathsf{T}},\dots,\mathbf{a}_{M}^{\mathsf{T}}]^{\mathsf{T}}$.
Subsequently, the rate of the $k$-th user decoding the $\ell$-th stream can be
written as

$\begin{split}&\underline{R}_{k,\ell}(\mathbf{a})=\log_{2}\Big{(}\varrho_{k}|(\mathbf{e}_{k,\ell}^{k,\ell})^{\mathsf{H}}\mathbf{a}|^{2}\\\
&+\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}(|(\mathbf{e}_{k,\ell}^{i,v})^{\mathsf{H}}\mathbf{a}|+\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2})^{2}+\sigma^{2}\Big{)}\\\
&-\log_{2}\Big{(}\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}(|(\mathbf{e}_{k,\ell}^{i,v})^{\mathsf{H}}\mathbf{a}|+\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2})^{2}+\sigma^{2}\Big{)},\end{split}$
(43)

where
$\mathbf{e}_{k,\ell}^{i,v}=\tilde{\mathbf{D}}_{i,v}^{\mathsf{H}}\hat{\mathbf{G}}^{\mathsf{H}}(\boldsymbol{\Theta})\tilde{\mathbf{c}}_{k,\ell}$.

Then, the gradient of $\underline{R}_{k,\ell}(\mathbf{a})$ with respect to
$\mathbf{a}$ is given in (44).

$\begin{split}\nabla\underline{R}_{k,\ell}(\mathbf{a})=\frac{2}{\ln
2}\bigg{[}&\frac{\varrho_{k}\mathbf{e}_{k,\ell}^{k,\ell}(\mathbf{e}_{k,\ell}^{k,\ell})^{\mathsf{H}}+\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}\big{(}1+\frac{\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2}}{|(\mathbf{e}_{k,\ell}^{k,\ell})^{\mathsf{H}}\mathbf{a}|}\big{)}\mathbf{e}_{k,\ell}^{i,v}(\mathbf{e}_{k,\ell}^{i,v})^{\mathsf{H}}}{\varrho_{k}|(\mathbf{e}_{k,\ell}^{k,\ell})^{\mathsf{H}}\mathbf{a}|^{2}+\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}(|(\mathbf{e}_{k,\ell}^{i,v})^{\mathsf{H}}\mathbf{a}|+\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2})^{2}+\sigma^{2}}\\\
&-\frac{\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}\big{(}1+\frac{\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2}}{|(\mathbf{e}_{k,\ell}^{k,\ell})^{\mathsf{H}}\mathbf{a}|}\big{)}\mathbf{e}_{k,\ell}^{i,v}(\mathbf{e}_{k,\ell}^{i,v})^{\mathsf{H}}}{\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}(|(\mathbf{e}_{k,\ell}^{i,v})^{\mathsf{H}}\mathbf{a}|+\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2})^{2}+\sigma^{2}}\bigg{]}\mathbf{a}.\end{split}$
(44)

Thus the gradient of (40a) with respect to $\mathbf{a}$ is given by

$\nabla_{\mathbf{a}}\underline{R}_{\text{sum}}(\mathbf{a})=\sum_{k=1}^{K}\sum_{\ell=1}^{L}\nabla_{\mathbf{a}}\underline{R}_{k,\ell}(\mathbf{a}).$
(45)

Similarly, we can obtain the gradient of $\mathcal{G}(\mathbf{a})$ with
respect to $\mathbf{a}$ as

$\begin{split}&\nabla_{\mathbf{a}}\mathcal{G}(\mathbf{a})=2\rho
K\Big{(}\sum_{k=1}^{K}\frac{1}{w_{k}^{2}}\sum_{\ell=1}^{L}\underline{R}_{k,\ell}(\mathbf{a})\sum_{m=1}^{L}\nabla_{\mathbf{a}}\underline{R}_{k,m}(\mathbf{a})\Big{)}\\\
&-2\Big{(}\sum_{k=1}^{K}\frac{1}{w_{k}}\sum_{\ell=1}^{L}\underline{R}_{k,\ell}(\mathbf{a})\Big{)}\Big{(}\sum_{i=1}^{K}\frac{1}{w_{i}}\sum_{m=1}^{L}\nabla_{\mathbf{a}}\underline{R}_{i,m}(\mathbf{a})\Big{)}.\end{split}$
(46)

Therefore, the gradient of the augmented Lagrangian function (41) with respect
to $\mathbf{a}$ can be written as

$\nabla_{\mathbf{a}}\mathcal{H}_{\gamma,\omega}(\mathbf{a})=\frac{\nabla_{\mathbf{a}}\underline{R}_{\text{sum}}(\mathbf{a})}{P_{\text{tot}}}-\gamma\nabla_{\mathbf{a}}\mathcal{G}(\mathbf{a})-\frac{1}{\omega}\
\mathcal{G}(\mathbf{a})\nabla_{\mathbf{a}}\mathcal{G}(\mathbf{a}).$ (47)

Hence, the update of $\mathbf{a}$ should follow

$\bar{\mathbf{a}}^{(t+1)}=\mathbf{a}^{(t)}+\alpha_{\mathbf{a}}^{(t)}\nabla_{\mathbf{a}}\mathcal{H}_{\gamma,\omega}(\mathbf{a}^{(t)}),$
(48)

and the projection of $\bar{\mathbf{a}}^{(t+1)}$ onto the constraint (16e) can
be performed as

$\mathbf{a}^{(t+1)}=\frac{1}{\sqrt{N_{\mathsf{T}}}}\text{diag}(|\bar{\mathbf{a}}^{(t+1)}|)^{-1}\bar{\mathbf{a}}^{(t+1)}.$
(49)

###  IV-C RIS Reflection Matrix

Next, assuming that the values of $\mathbf{D},\mathbf{A}$ and
$\tilde{\mathbf{C}}$ are fixed, the RIS reflection matrix optimization sub-
problem can be stated as

$\displaystyle\max_{\boldsymbol{\Theta}}\ \ \ $
$\displaystyle\underline{\eta}(\boldsymbol{\Theta}),$ (50a) s.t.
$\displaystyle\eqref{F_const},\ \eqref{ps_con_EEf},$ (50b)

with the following augmented Lagrangian function:

$\mathcal{H}_{\gamma,\omega}(\boldsymbol{\Theta})=\frac{1}{P_{\text{tot}}}\underline{R}_{\text{sum}}(\boldsymbol{\Theta})-\gamma\mathcal{G}(\boldsymbol{\Theta})-\frac{1}{2\omega}\mathcal{G}^{2}(\boldsymbol{\Theta}).$
(51)

To find the gradient of $\mathcal{H}_{\gamma,\omega}(\boldsymbol{\Theta})$
with respect to $\boldsymbol{\Theta}$, we use the fact that the vectorization
of a scalar term is the term itself to re-express the terms
$\tilde{\mathbf{c}}_{k,\ell}^{\mathsf{H}}\hat{\mathbf{G}}_{k}(\boldsymbol{\Theta})\mathbf{A}\mathbf{d}_{i,v}$
as

$\tilde{\mathbf{c}}_{k,\ell}^{\mathsf{H}}\hat{\mathbf{G}}_{k}(\boldsymbol{\Theta})\mathbf{A}\mathbf{d}_{i,v}=\underbrace{\big{(}(\mathbf{A}\mathbf{d}_{i,v})^{\mathsf{T}}\otimes\tilde{\mathbf{c}}_{k,\ell}^{\mathsf{H}}\big{)}\hat{\mathbf{H}}_{k}}_{(\mathbf{f}_{k,\ell}^{i,v})^{\mathsf{H}}}\boldsymbol{\theta},$
(52)

which can be obtained by applying the identity
$\text{vec}(\boldsymbol{\Omega}_{1}\boldsymbol{\Omega}_{2}\boldsymbol{\Omega}_{3})=(\boldsymbol{\Omega}_{3}^{\mathsf{T}}\otimes\boldsymbol{\Omega}_{1})\text{vec}(\boldsymbol{\Omega}_{2})$,
which holds for any general matrices $\boldsymbol{\Omega}_{1}$,
$\boldsymbol{\Omega}_{2}$ and $\boldsymbol{\Omega}_{3}$ with appropriate
dimensions [30].

Thus, we can write the rate of the $k$-th user decoding the $\ell$-th stream
as

$\begin{split}&\underline{R}_{k,\ell}(\boldsymbol{\theta})=\log_{2}\Big{(}\varrho_{k}|(\mathbf{f}_{k,\ell}^{k,\ell})^{\mathsf{H}}\boldsymbol{\theta}|^{2}\\\
&+\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}\big{(}|(\mathbf{f}_{k,\ell}^{i,v})^{\mathsf{H}}\boldsymbol{\theta}|+\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2}\big{)}^{2}+\sigma^{2}\Big{)}\\\
&-\log_{2}\Big{(}\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}\big{(}|(\mathbf{f}_{k,\ell}^{i,v})^{\mathsf{H}}\boldsymbol{\theta}|+\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2}\big{)}^{2}+\sigma^{2}\Big{)}.\end{split}$
(53)

Then, the gradient of $\underline{R}_{k,\ell}(\boldsymbol{\theta})$ with
respect to $\boldsymbol{\theta}$ is given by (54), shown below,

$\begin{split}\nabla_{\boldsymbol{\theta}}\underline{R}_{k,\ell}(\boldsymbol{\theta})=\frac{2}{\ln
2}\bigg{[}&\frac{\varrho_{k}\mathbf{f}_{k,\ell}^{k,\ell}(\mathbf{f}_{k,\ell}^{k,\ell})^{\mathsf{H}}+\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}\big{(}1+\frac{\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2}}{|(\mathbf{f}_{k,\ell}^{i,v})^{\mathsf{H}}\boldsymbol{\theta}|}\big{)}\mathbf{f}_{k,\ell}^{i,v}(\mathbf{f}_{k,\ell}^{i,v})^{\mathsf{H}}}{\varrho_{k}|(\mathbf{f}_{k,\ell}^{k,\ell})^{\mathsf{H}}\boldsymbol{\theta}|^{2}+\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}\big{(}|(\mathbf{f}_{k,\ell}^{i,v})^{\mathsf{H}}\boldsymbol{\theta}|+\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2}\big{)}^{2}+\sigma^{2}}\\\
&-\frac{\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}\big{(}1+\frac{\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2}}{|(\mathbf{f}_{k,\ell}^{i,v})^{\mathsf{H}}\boldsymbol{\theta}|}\big{)}\mathbf{f}_{k,\ell}^{i,v}(\mathbf{f}_{k,\ell}^{i,v})^{\mathsf{H}}}{\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{v=1}^{L}}\big{(}|(\mathbf{f}_{k,\ell}^{i,v})^{\mathsf{H}}\boldsymbol{\theta}|+\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2}\big{)}^{2}+\sigma^{2}}\bigg{]}\boldsymbol{\theta}.\end{split}$
(54)

and that of $\underline{R}_{\text{sum}}(\boldsymbol{\theta})$ is

$\nabla_{\boldsymbol{\theta}}\underline{R}_{\text{sum}}(\boldsymbol{\theta})=\sum_{k=1}^{K}\sum_{\ell=1}^{L}\nabla_{\boldsymbol{\theta}}\underline{R}_{k,\ell}(\boldsymbol{\theta}).$
(55)

Furthermore, the gradient of $\mathcal{G}(\boldsymbol{\theta})$ with respect
to $\boldsymbol{\theta}$ is obtained as

$\begin{split}&\nabla_{\boldsymbol{\theta}}\mathcal{G}({\boldsymbol{\theta}})=2\rho
K\Big{(}\sum_{k=1}^{K}\frac{1}{w_{k}^{2}}\sum_{\ell=1}^{L}\underline{R}_{k,\ell}({\boldsymbol{\theta}})\sum_{m=1}^{L}\nabla_{\boldsymbol{\theta}}\underline{R}_{k,m}({\boldsymbol{\theta}})\Big{)}\\\
&-2\Big{(}\sum_{k=1}^{K}\frac{1}{w_{k}}\sum_{\ell=1}^{L}\underline{R}_{k,\ell}({\boldsymbol{\theta}})\Big{)}\Big{(}\sum_{i=1}^{K}\frac{1}{w_{i}}\sum_{m=1}^{L}\nabla_{\boldsymbol{\theta}}\underline{R}_{i,m}({\boldsymbol{\theta}})\Big{)}.\end{split}$
(56)

Therefore, the gradient $\mathcal{H}_{\gamma,\omega}$ with respect to
$\boldsymbol{\theta}$ is

$\nabla_{\boldsymbol{\theta}}\mathcal{H}_{\gamma,\omega}(\boldsymbol{\theta})=\frac{\nabla_{\boldsymbol{\theta}}\underline{R}_{\text{sum}}(\boldsymbol{\theta})}{P_{\text{tot}}}-\gamma\nabla_{\boldsymbol{\theta}}\mathcal{G}(\boldsymbol{\theta})-\frac{1}{\omega}\mathcal{G}(\boldsymbol{\theta})\nabla_{\boldsymbol{\theta}}\mathcal{G}(\boldsymbol{\theta}).$
(57)

Therefore, the update of $\boldsymbol{\theta}$ should follow

$\bar{\boldsymbol{\theta}}^{(t+1)}=\boldsymbol{\theta}^{(t)}+\alpha_{\boldsymbol{\theta}}^{(t)}\nabla_{\boldsymbol{\theta}}\mathcal{H}_{\gamma,\omega}(\boldsymbol{\theta}^{(t)}),$
(58)

and the projection onto the constraint (16d) can be performed as

$\boldsymbol{\theta}^{(t+1)}=\text{diag}(|\bar{\boldsymbol{\theta}}^{(t+1)}|)^{-1}\bar{\boldsymbol{\theta}}^{(t+1)}.$
(59)

###  IV-D Receive Combiners

Assuming that the values of $\mathbf{D},\mathbf{A}$ and $\boldsymbol{\Theta}$
are fixed, the combiner optimization sub-problem can be written as

$\displaystyle\max_{\tilde{\mathbf{C}}}\ \ \ $
$\displaystyle\underline{\eta}(\tilde{\mathbf{C}}),$ (60a) s.t.
$\displaystyle\eqref{F_const},\ \eqref{lb_constraint},$ (60b)

with the following augmented Lagrangian function:

$\mathcal{H}_{\gamma,\omega}(\tilde{\mathbf{C}})=\frac{1}{P_{\text{tot}}}\underline{R}_{\text{sum}}(\tilde{\mathbf{C}})-\gamma\mathcal{G}(\tilde{\mathbf{C}})-\frac{1}{2\omega}\
\mathcal{G}^{2}(\tilde{\mathbf{C}}).$ (61)

We start by expressing $\underline{R}_{k,\ell}(\tilde{\mathbf{c}}_{k,\ell})$
as

$\begin{split}&\underline{R}_{k,\ell}(\tilde{\mathbf{c}}_{k,\ell})=\log_{2}\bigg{(}\varrho_{k}|(\mathbf{J}_{k}^{k,\ell})^{\mathsf{H}}\tilde{\mathbf{c}}_{k,\ell}|^{2}\\\
&+\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{{v=1}}^{L}}\big{(}|(\mathbf{J}_{k}^{i,v})^{\mathsf{H}}\tilde{\mathbf{c}}_{k,\ell}|+\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2}\big{)}^{2}+\sigma^{2}\bigg{)}\\\
&-\log_{2}\bigg{(}\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{{v=1}}^{L}}\big{(}|(\mathbf{J}_{k}^{i,v})^{\mathsf{H}}\tilde{\mathbf{c}}_{k,\ell}|+\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2}\big{)}^{2}\bigg{)},\end{split}$
(62)

with
$\mathbf{J}_{k}^{i,v}=\hat{\mathbf{G}}_{k}(\boldsymbol{\Theta})\mathbf{A}\mathbf{d}_{i,v}$.

The gradient of $\underline{R}_{k,\ell}(\tilde{\mathbf{c}}_{k,\ell})$ with
respect to $\tilde{\mathbf{c}}_{k,\ell}$can then be obtained as expressed in
(63) on the next page,

$\begin{split}\nabla_{\tilde{\mathbf{c}}_{k,\ell}}\underline{R}_{k,\ell}(\tilde{\mathbf{c}}_{k,\ell})=\frac{2}{\ln
2}\bigg{[}&\frac{\varrho_{k}\mathbf{J}_{k}^{i,v}(\mathbf{J}_{k}^{i,v})^{\mathsf{H}}+\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{{v=1}}^{L}}\big{(}1+\frac{\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2}}{|(\mathbf{J}_{k}^{i,v})^{\mathsf{H}}\tilde{\mathbf{c}}_{k,\ell}}\big{)}\mathbf{J}_{k}^{i,v}(\mathbf{J}_{k}^{i,v})^{\mathsf{H}}}{\varrho_{k}|(\mathbf{J}_{k}^{k,\ell})^{\mathsf{H}}\tilde{\mathbf{c}}_{k,\ell}|^{2}+\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{{v=1}}^{L}}\big{(}|(\mathbf{J}_{k}^{i,v})^{\mathsf{H}}\tilde{\mathbf{c}}_{k,\ell}|+\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2}\big{)}^{2}+\sigma^{2}}\\\
&-\frac{\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{{v=1}}^{L}}\big{(}1+\frac{\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2}}{|(\mathbf{J}_{k}^{i,v})^{\mathsf{H}}\tilde{\mathbf{c}}_{k,\ell}}\big{)}\mathbf{J}_{k}^{i,v}(\mathbf{J}_{k}^{i,v})^{\mathsf{H}}}{\underset{[i,v]\neq[k,\ell]}{\sum_{i=1}^{K}\sum_{{v=1}}^{L}}\big{(}|(\mathbf{J}_{k}^{i,v})^{\mathsf{H}}\tilde{\mathbf{c}}_{k,\ell}|+\delta_{k}\sqrt{N_{\text{RIS}}}||\mathbf{d}_{i,v}||_{2}\big{)}^{2}+\sigma^{2}}\bigg{]}\tilde{\mathbf{c}}_{k,\ell},\end{split}$
(63)

and that of $\underline{R}_{\text{sum}}(\tilde{\mathbf{C}})$ with respect to
$\tilde{\mathbf{C}}$ is

$\nabla_{\tilde{\mathbf{C}}}\underline{R}_{\text{sum}}(\tilde{\mathbf{C}})=[\nabla_{\tilde{\mathbf{c}}_{k,\ell}}\underline{R}_{1,1}(\tilde{\mathbf{c}}_{1,1}),\dots,\nabla_{\tilde{\mathbf{c}}_{k,\ell}}\underline{R}_{K,L}(\tilde{\mathbf{c}}_{K,L})].\vspace{-0.5cm}$
(64)

Additionally, the gradient of $\mathcal{G}(\tilde{\mathbf{C}})$ with respect
to $\tilde{\mathbf{C}}$ is found as

$\begin{split}&\nabla_{\tilde{\mathbf{C}}}\mathcal{G}({\tilde{\mathbf{C}}})=2\rho
K\Big{(}\sum_{k=1}^{K}\frac{1}{w_{k}^{2}}\sum_{\ell=1}^{L}\underline{R}_{k,\ell}({\tilde{\mathbf{c}}}_{k,\ell})\\\
&\times\sum_{m=1}^{L}\nabla_{\mathbf{D}}\underline{R}_{k,m}({\tilde{\mathbf{c}}}_{k,m})\Big{)}\\\
&-2\Big{(}\sum_{k=1}^{K}\frac{1}{w_{k}}\sum_{\ell=1}^{L}\underline{R}_{k,\ell}({\tilde{\mathbf{c}}}_{k,\ell})\Big{)}\Big{(}\sum_{i=1}^{K}\frac{1}{w_{i}}\sum_{m=1}^{L}\nabla_{\mathbf{D}}\underline{R}_{i,m}({\tilde{\mathbf{c}}}_{i,m})\Big{)}.\end{split}$
(65)

Thus, the gradient of $\mathcal{H}_{\gamma,\omega}(\tilde{\mathbf{C}})$ with
respect to $\tilde{\mathbf{C}}$ is

$\nabla_{\tilde{\mathbf{C}}}\mathcal{H}_{\gamma,\omega}(\tilde{\mathbf{C}})=\frac{\nabla_{\tilde{\mathbf{C}}}\underline{R}_{\text{sum}}(\tilde{\mathbf{C}})}{P_{\text{tot}}}-\gamma\nabla_{\tilde{\mathbf{C}}}\mathcal{G}(\tilde{\mathbf{C}})-\frac{1}{\omega}\mathcal{G}(\tilde{\mathbf{C}})\nabla_{\tilde{\mathbf{C}}}\mathcal{G}(\tilde{\mathbf{C}}).$
(66)

Thus, the update of $\tilde{\mathbf{C}}$ should follow

$\bar{\mathbf{C}}^{(t+1)}=\tilde{\mathbf{C}}^{(t)}+\alpha_{\mathbf{C}}^{(t)}\nabla_{\tilde{\mathbf{C}}}\mathcal{H}_{\gamma,\omega}(\tilde{\mathbf{C}}^{(t)}),$
(67)

and the projection onto the constraint (16f) can be done by normalizing each
column of $\tilde{\mathbf{C}}$ to have unit norm, i.e.,

$\tilde{\mathbf{c}}_{k,\ell}^{(t+1)}=\bar{\mathbf{c}}_{k,\ell}^{(t+1)}/||\bar{\mathbf{c}}_{k,\ell}^{(t+1)}||_{2},\
\forall\ell,k.$ (68)

###  IV-E The Slack Variable $\mu$

The optimization problem with respect to the slack variable $\mu$ is

$\displaystyle\max_{\mu\geq 0}\ \ \ $
$\displaystyle\mathcal{H}_{\gamma,\omega}(\mu).$ (69a)

To find the optimal value of $\mu$, we find the partial derivative
$\frac{\partial\mathcal{H}_{\gamma,\omega}(\mu)}{\partial\mu}$ and equate it
to zero to obtain

$\frac{\partial\mathcal{H}_{\gamma,\omega}(\mu)}{\partial\mu}=-\gamma-\frac{1}{\omega}\mathcal{G}(\mu)=0,$
(70)

which has a solution of

$\begin{split}\bar{\mu}&=\Big{(}\sum_{k=1}^{K}\frac{1}{w_{k}}\sum_{\ell=1}^{L}\underline{R}_{k,\ell}\Big{)}^{2}-\rho
K\sum_{k=1}^{K}\Big{(}\frac{1}{w_{k}}\sum_{\ell=1}^{L}\underline{R}_{k,\ell}\Big{)}^{2}-\gamma\omega.\end{split}$
(71)

Then the update of $\mu$ after projecting onto the feasible set is

$\mu^{(t+1)}=\max\big{(}0,\bar{\mu}\big{)}.$ (72)

Algorithm 1 summarizes the proposed penalty-based projected gradient ascent
approach.

Algorithm 1 AO algorithm based on projected gradient ascent and PDD method

1:Input: $\rho\in[0,1]$, $\gamma\geq 0$, $\omega>0$, $\psi\in[0,1]$,
$\mathbf{D}^{(0)}$, $\mathbf{A}^{(0)}$, $\boldsymbol{\Theta}^{(0)}$,
$\tilde{\mathbf{C}}^{(0)}$, $\mu^{(0)}=0$, $\underline{\eta}^{*}$,
$\rho\in[0,1]$, $\gamma$, $\omega$, $t=0$, and $\varepsilon>0$.

2:While $\mu^{(t)}>0$ do {

3:do {

4:Update
$\bar{\mathbf{D}}^{(t+1)}\leftarrow\mathbf{D}^{(t)}+\alpha_{\mathbf{D}}^{(t)}\nabla_{\mathbf{D}}\mathcal{H}_{\gamma,\omega}(\mathbf{D}^{(t)})$
using (38).

5:Project onto the constraint (16c) using (39).

6:

7:Update
$\bar{\mathbf{a}}^{(t+1)}\leftarrow\mathbf{a}^{(t)}+\alpha_{\mathbf{a}}^{(t)}\nabla_{\mathbf{a}}\mathcal{H}_{\gamma,\omega}(\mathbf{a}^{(t)})$
using (48).

8:Project onto the constraint (16e) using (49).

9:

10:Update
$\bar{\boldsymbol{\theta}}^{(t+1)}\leftarrow\boldsymbol{\theta}^{(t)}+\alpha_{\boldsymbol{\theta}}^{(t)}\nabla_{\boldsymbol{\theta}}\mathcal{H}_{\gamma,\omega}(\boldsymbol{\theta}^{(t)})$
using (58).

11:Project onto the constraint (16d) using (59).

12:

13:Update
$\bar{\mathbf{C}}^{(t+1)}\leftarrow\tilde{\mathbf{C}}^{(t)}+\alpha_{\mathbf{C}}^{(t)}\nabla\mathcal{H}_{\gamma,\omega}(\tilde{\mathbf{C}}^{(t)})$
using (67).

14:Project onto the constraint (16f) using (68).

15:

16:Update $\mu^{(t+1)}$ using (72).

17:

18:Update $t\leftarrow t+1.$

19: }

20:until | $\mathcal{H}_{\gamma,\omega}(\mathbf{D}^{(t)},\mathbf{A}^{(t)},\boldsymbol{\Theta}^{(t)},\tilde{\mathbf{C}}^{(t)},\mu^{(t)})-\mathcal{H}_{\gamma,\omega}($ $\mathbf{D}^{(t-1)},\mathbf{A}^{(t-1)},\boldsymbol{\Theta}^{(t-1)},\tilde{\mathbf{C}}^{(t-1)},\mu^{(t-1)})\big{|}\leq\varepsilon$.

21:Update
$\gamma\leftarrow\gamma+\frac{1}{\omega}\mathcal{H}_{\gamma,\omega}(\mathbf{D}^{(t)},\mathbf{A}^{(t)},\boldsymbol{\Theta}^{(t)},\tilde{\mathbf{C}}^{(t)},$
$\mu^{(t)})$.

22:Update $\omega\leftarrow\psi\omega$.

23:}

24:Output:
$\underline{\eta}^{*}=\underline{\eta}(\mathbf{D}^{(t)},\mathbf{A}^{(t)},\boldsymbol{\Theta}^{(t)},\tilde{\mathbf{C}}^{(t)})$.

###  IV-F Convergence

The employed algorithms (i.e., AO, projected gradient ascent, and the PDD
method) are proven to converge to a locally optimal solution, as was
demonstrated in [31, 32, 29]. Due to the nature of the targeted optimization
problem, finding a global solution poses a challenge. However, a satisfactory
locally optimal solution can be obtained by solving the optimization problem
multiple times with different random initializations, as will be illustrated
in Section V-A.

###  IV-G Complexity Analysis

In this subsection, we provide an analysis of the computational complexity of
the proposed method, defined as the number of required complex-valued
multiplications. Each variable update and projection for Algorithm 1 is
dominated by specific vector multiplications, including inner and outer
products. Table I outlines the complexity of updating and projecting each
variable.

Table I: Computational complexity. Variable |  Update complexity |  Projection complexity
---|---|---
$\mathbf{D}$ |  $\mathcal{O}\big{(}KLM(MLN_{\mathsf{T}}N_{\text{RIS}}+N_{\text{RIS}}^{2}+N_{\text{RIS}}N_{\mathsf{R}}+KL)\big{)}$ |  $\mathcal{O}\big{(}K^{2}L^{2}M\big{)}$
$\mathbf{A}$ |  $\mathcal{O}\big{(}K^{2}L^{2}MN_{\mathsf{T}}(MN_{\mathsf{T}}N_{\text{RIS}}+N_{\text{RIS}}^{2}+N_{\text{RIS}}N_{\mathsf{R}})\big{)}$ |  $\mathcal{O}\big{(}MN_{\mathsf{T}}\big{)}$
$\boldsymbol{\Theta}$ |  $\mathcal{O}\big{(}K^{2}L^{2}N_{\text{RIS}}(N_{\mathsf{R}}+M^{2}N_{\mathsf{T}}+N_{\text{RIS}}\big{)}$ |  $\mathcal{O}\big{(}N_{\text{RIS}}\big{)}$
$\tilde{\mathbf{C}}$ |  $\mathcal{O}\big{(}K^{2}L^{2}N_{\mathsf{R}}(N_{\text{RIS}}^{2}+M^{2}N_{\mathsf{T}}+N_{\mathsf{R}})\big{)}$ |  $\mathcal{O}\big{(}KLN_{\mathsf{R}}\big{)}$

In practical systems, the number of RIS elements is significantly larger than
the number the number of BS antennas, the number of UE antennas, and the
number of data streams i.e.,
$N_{\text{RIS}}\gg\max(MN_{\mathsf{T}},N_{\mathsf{R}},LK)$. In this scenario,
the computational complexity of one inner iteration of Algorithm 1 is
primarily dominated by the term
$\mathcal{O}\big{(}L^{2}K^{2}N_{\text{RIS}}^{2}(MN_{\mathsf{T}}+N_{\mathsf{R}})\big{)}$.

##  V Numerical Simulations

In this section, we conduct extensive simulations to demonstrate the
performance achievable via the proposed method. We consider two state-of-the-
art benchmark schemes for comparison, the EE-only maximization scheme of [12]
(both without and with rate constraints), and the EE fairness maximization
scheme of [23].

In our simulations, we average the performance metric over 1000 independent
channel realizations, with small-scale fading, user locations, and QoS weights
randomly varied within confined ranges. Table II provides a summary of the
simulation parameters; these are adopted in all simulations unless stated
otherwise.

Table II: Simulation parameters. Parameter |  Value |  Parameter |  Value
---|---|---|---
$M$ |  8 |  $P_{\theta}$ |  1 dBm
$N_{\mathsf{T}}$ |  4 |  $\xi$ |  1.2
$N_{\text{RIS}}$ |  64 |  $\rho$ |  0.75
$K$ |  4 |  $\varepsilon$ |  $1\times 10^{-3}$
$L$ |  2 |  $\sigma^{2}$ |  -37 dBm
$N_{\mathsf{R}}$ |  4 |  $w_{k}$ range |  $[1,5]$
Carrier frequency |  $28\,\mathrm{GHz}$ |  Location of the reference AE at the BS |  $(0,0,10)\,\mathrm{m}$
Total bandwidth |  $200\,\mathrm{MHz}$ |  Location of the center of the RIS |  $(15,-15,5)\,\mathrm{m}$
$P_{\text{max}}$ |  40 dBm |  $k$-th UE’s location range |  $[(15,-15,0)\,\mathrm{m},$ $(30,15,2)\,\mathrm{m}]$
$P_{\text{BS}}$ |  9 dBW |  $k$-th user’s imperfect CSI bound |  $0.2||\hat{\mathbf{H}}_{k}||_{\mathsf{F}}$
$P_{\text{UE},k}$ |  5 dBm |  Rate constraints for [12] |  $R_{k}\geq 0.1w_{k}\log_{2}\big{(}\frac{P_{\text{max}}}{K\sigma^{2}}\big{)}$
$P_{\text{RF,T}}$ |  1 dBm |  |

###  V-A Convergence

|
---|---
(a) Perfect CSI. | (b) Imperfect CSI.
Figure 2: EE (solid) and the augmented Lagrangian function (dashed) versus
iteration number for the cases of perfect and imperfect CSI.

We first examine the convergence of the algorithm by plotting the EE and the
augmented Lagrangian function versus the iteration number for both perfect and
imperfect CSI cases in Fig. 2. We use an initial step size of one to update
each variable and progressively decrease it by half in the case of
overshooting (i.e., when we observe a decrease in the objective function). We
initiate the process by setting the values of $\gamma$ and $\omega$ to 0 and
10, respectively. Using eight random initial points that satisfy the
constraint set, we observe convergence towards two sub-optimal points due to
the non-convexity of the problem. In this initial stage, the weight of the
constraint in the augmented Lagrangian function is small. Therefore, the
algorithm primarily focuses on optimizing the EE part of the objective
function. Since the value of the EE in this case is almost equal to that of
the augmented Lagrangian function, we omit the plot of the augmented
Lagrangian function. Subsequently, the parameter $\omega$ is reduced by a
factor of 10 to accentuate the constraint more in the augmented Lagrangian
function until it is satisfied.

At the final iteration, the augmented Lagrangian function converges to the
value of the EE objective function. This convergence is expected, as both
functions should have the same value when the constraint (16b) is satisfied.
The algorithm eventually converges to two EE values (179 and 145
Mbit/sec/Joule for the perfect CSI scenario, and 149 and 128 Mbit/sec/Joule
for the imperfect CSI scenario). To ensure convergence to a satisfactory
solution, we utilize the random initialization technique in all simulations by
employing five initial points and selecting the highest final EE value.

###  V-B Effect of Varying the Fairness Requirement

Figure 3: Average EE versus Jain’s fairness index for different values of
$\rho$.

Fig. 3 demonstrates the EE-fairness trade-off achieved by varying the value of
$\rho$. The general trend indicates that as $\rho$ increases, the EE
decreases, and the Jain’s fairness index increases. This relationship arises
because $\rho$ controls the fairness requirement of the system at the cost of
reducing the EE value. At one extreme, perfect fairness (Jain’s fairness index
of one) can be achieved when $\rho=1$, but this comes at the expense of
reducing the EE by 63% of its optimal value when no fairness is imposed. At
the other extreme, the maximum EE can be achieved with very poor fairness. The
value of $\rho$ can be tuned to achieve a favorable and flexible trade-off
between these two extremes.

Compared to the existing EE maximization method of [12], the proposed
algorithm can achieve a better EE value with the same fairness index. This
improvement is attributed to our proposed system model and transmission
scheme, which are more suitable for EE applications than that of [12]. While
the poor fairness in [12] can be enhanced by incorporating QoS constraints, it
comes with a significant reduction in the EE. On the other hand, the fairness
algorithm proposed in [23] can achieve a point close to the Pareto-optimal
front generated by varying $\rho$ in our proposed algorithm. Nevertheless, the
proposed method offers the flexibility of tuning $\rho$ to achieve the desired
operating point of EE and fairness based on available resources and specified
QoS requirements.

###  V-C Effect of Varying the CSI Error Bound

Figure 4: Behavior of the average EE with increasing CSI error bound.

Fig. 4 shows the effect of increasing the CSI error bound, where the CSI error
bound of each user is modeled as

$\delta_{k}=\beta||\hat{\mathbf{H}}_{k}||_{\mathsf{F}},$ (73)

for some $\beta\in[0,0.2]$. We plot three curves: the true EE given in (13) by
optimizing the lower bound in (15), the lower bound on the EE given in (15),
and the EE in the case when the BS assumes no error exists in the CSI.

It can be noticed that, compared to the case where the BS assumes perfect CSI
knowledge, optimizing the lower bound instead of the original function can
offer improvements in EE. However, this improvement comes at an additional
cost in complexity, as more terms must be included in the gradient
calculations. In conclusion, if the imperfect CSI bound is small, it might be
more efficient for the BS to assume perfect CSI knowledge, as the improvement
is not significant. Conversely, optimizing the EE lower bound offers a notable
improvement for scenarios where the CSI uncertainty is relatively high.

###  V-D Effect of Varying the Transmit Power Budget

|
---|---
(a) Average EE. | (b) Average Jain’s fairness index.
Figure 5: Average EE and Jain’s fairness index as the transmit power budget
increases.

Fig. 5 shows the average EE and Jain’s fairness index as the maximum transmit
power budget increases. It can be noted that both EE and fairness saturate
when the transmit power budget exceeds 40 dBm. This indicates that the extra
power is not being used as it does not help in optimizing the EE or fairness.
This can be attributed to the fact that the numerator of the EE objective
function is a logarithmic function of the transmit power while the denominator
is a linear function of the transmit power. In other terms, the rate of
increase of the numerator (sum-rate) decreases with increasing
$P_{\text{max}}$, while that of the denominator (total power consumption) is
increasing in a constant rate.

Furthermore, it can be noted that $\rho=0$ in the proposed method can achieve
a better EE than the state-of-art methods for three main reasons. First, the
AoSA architecture is more energy-efficient than the fully-connected one. This
conclusion matches observation made in [19], where the AoSA and fully-
connected architectures are compared from this perspective. Second, the
proposed system has a greater capability of increasing the EE as it features
two additional sets of variables: the analog precoder at the BS and the
combiners at the UEs. Third, unlike the benchmark methods, ZF and orthogonal
transmission have not been assumed. While orthogonal transmission clearly
limits the system capabilities, ZF is not ideal in the case of correlated,
ill-conditioned, and rank-deficient mmWave channels. Nevertheless, considering
only EE leads to very poor fairness, as in most practical cases all resources
are allocated to serve one user only. This issue can be resolved by targeting
fairness optimization directly as in [23], however, the resulting EE is poorer
than that achieved by other methods. The proposed approach can achieve a
better EE and fairness trade-off than these existing approaches, as optimizing
the EE with a Jain’s fairness index requirement of 0.75 can lead to highly
favorable values for both the EE and user fairness. On the one hand, although
incorporating rate constraints can help to improve fairness, it can degrade
the system EE significantly when these constraints depend on the transmit
power budget. Nevertheless, setting the values for the rate constraints
independently of the transmit power budget may lead to feasibility issues for
the optimization problem.

###  V-E Effect of Varying the Number of RIS elements

Figure 6: Average EE as the number of RIS elements increases.

In Fig. 6, we plot the average EE as the number of RIS elements increases for
two values of the power consumption at each RIS element:
$P_{\theta}=1\,\mathrm{dBm}$ and $P_{\theta}=10\,\mathrm{dBm}$. We can notice
that adding extra RIS elements provides a substantial increase in the EE for
all cases when the number of RIS elements is small. However, as the number of
RIS elements becomes larger than 64, this increase begins to slow down in the
case of $P_{\theta}=1\,\mathrm{dBm}$, while the average EE starts to decrease
for $P_{\theta}=10\,\mathrm{dBm}$. This eventual decrease in the EE occurs due
to the fact that although more RIS elements offer more degrees of freedom to
enhance the system performance, it provides a logarithmic increase in the
numerator of the EE but a linear increase in the denominator.

##  VI Conclusion

This paper has presented a method for robust design of the hybrid analog and
digital precoder at the BS, RIS reflection matrix, and digital combiners at
the UEs in order to maximize the system EE while imposinga specified Jain’s
fairness index in RIS-assisted mmWave systems with imperfect channel state
information. To achieve this, a lower bound based on the triangle and Cauchy-
Schwarz inequalities was employed. The penalty dual decomposition method was
then used to address the challenging fairness constraint, and the projected
gradient ascent method was utilized to obtain the solution to the optimization
problem. The simulation results show that the proposed method can provide a
very good tradeoff between EE and fairness, while also offering a high degree
of flexibility in tuning the EE and user fairness in order to prioritize one
of these metrics over the other.

##  Appendix A Proof of Proposition 1

The signal-to-interference-plus-noise ratio (SINR) in (11) consists of three
components, the desired signal term, the interference term, and the noise
term.

To obtain a lower bound on the desired signal term, we seek the error matrix
that minimizes that term by solving the following convex optimization problem:

$\displaystyle\min_{\boldsymbol{\Lambda}_{k}}\ \ $
$\displaystyle\mathcal{D}_{k}(\boldsymbol{\Lambda}_{k})=\big{|}\mathbf{c}_{k,\ell}^{\mathsf{H}}\big{(}\hat{\mathbf{G}}_{k}(\boldsymbol{\Theta})+\boldsymbol{\Lambda}_{k}\big{)}\mathbf{A}\mathbf{d}_{k,\ell}\big{|}^{2},$
(74a) s.t.
$\displaystyle||\boldsymbol{\Delta}_{k}||_{\text{F}}\leq\delta_{k}.$ (74b)

We can rewrite the objective function (74a) as

$\begin{split}\mathcal{D}_{k}(\boldsymbol{\Delta}_{k})=\big{|}\big{(}(\mathbf{A}\mathbf{d}_{k,\ell})^{\mathsf{T}}\otimes\mathbf{c}_{k,\ell}^{\mathsf{H}}\big{)}\big{(}\hat{\mathbf{H}}_{k}+\boldsymbol{\Delta}_{k}\big{)}\boldsymbol{\theta}\big{|}^{2}.\end{split}$
(75)

We observe that (75) attains its minimum when
$\boldsymbol{\Delta}_{k}=-\hat{\mathbf{H}}_{k}$. By projecting this solution
onto the feasible set, the optimal solution is found as
$\boldsymbol{\Delta}_{k}^{*}=-\min\big{(}1,\delta_{k}/||\hat{\mathbf{H}}_{k}||_{\mathsf{F}}\big{)}\hat{\mathbf{H}}_{k}$.
Assuming the employed channel estimation algorithm is reliable enough such
that $||\hat{\mathbf{H}}_{k}||_{\mathsf{F}}\geq\delta_{k}$, the optimal
solution becomes
$\boldsymbol{\Delta}_{k}^{*}=-\delta_{k}\hat{\mathbf{H}}_{k}/||\hat{\mathbf{H}}_{k}||_{\mathsf{F}}$.
Substituting this solution into the objective function (74a), we obtain a
lower bound on the desired signal component as:

$\mathcal{D}_{k}(\boldsymbol{\Delta}_{k})\geq\Big{(}1-\frac{\delta_{k}}{||\hat{\mathbf{H}}_{k}||_{\mathsf{F}}}\Big{)}^{2}\big{|}\mathbf{c}_{k,\ell}^{\mathsf{H}}\hat{\mathbf{G}}_{k}(\boldsymbol{\Theta})\mathbf{A}\mathbf{d}_{k,\ell}\big{|}^{2}.$
(76)

For the interference term, we use the triangle inequality followed by the
Cauchy–Schwarz inequality to obtain

$\begin{split}&\big{|}\mathbf{c}_{k,\ell}^{\mathsf{H}}\big{(}\hat{\mathbf{G}}_{k}(\boldsymbol{\Theta})+\boldsymbol{\Lambda}_{k}\big{)}\mathbf{A}\mathbf{d}_{i,v}\big{|}\\\
&\leq|\mathbf{c}_{k,\ell}^{\mathsf{H}}\hat{\mathbf{G}}_{k}(\boldsymbol{\Theta})\mathbf{A}\mathbf{d}_{i,v}|+|\mathbf{c}_{k,\ell}^{\mathsf{H}}\boldsymbol{\Lambda}_{k}\mathbf{A}\mathbf{d}_{i,v}|\\\
&\leq|\mathbf{c}_{k,\ell}^{\mathsf{H}}\hat{\mathbf{G}}_{k}(\boldsymbol{\Theta})\mathbf{A}\mathbf{d}_{i,v}|+||\mathbf{c}_{k,\ell}||_{2}||\boldsymbol{\Lambda}_{k}||_{\mathsf{F}}||\mathbf{A}\mathbf{d}_{i,v}||_{2}.\end{split}$
(77)

Using the facts
$||\boldsymbol{\Lambda}_{k}||_{\mathsf{F}}\leq||\boldsymbol{\Delta}_{k}||_{\mathsf{F}}||\boldsymbol{\theta}_{k}||_{2}\leq\delta_{k}\sqrt{N_{\mathsf{RIS}}}$
and $||\mathbf{A}\mathbf{d}_{i,v}||_{2}=||\mathbf{d}_{i,v}||_{2}$, we can
obtain

$\begin{split}|\mathbf{c}_{k,\ell}^{\mathsf{H}}\mathbf{G}_{k}(\boldsymbol{\Theta})\mathbf{A}\mathbf{d}_{i,v}|\leq&|\mathbf{c}_{k,\ell}^{\mathsf{H}}\hat{\mathbf{G}}_{k}(\boldsymbol{\Theta})\mathbf{A}\mathbf{d}_{i,v}|\\\
&+\delta_{k}\sqrt{N_{\mathsf{RIS}}}||\mathbf{c}_{k,\ell}||_{2}||\mathbf{d}_{i,v}||_{2}.\end{split}$
(78)

The noise term can be upper bounded as

$\mathbb{E}\\{|\mathbf{c}_{k,\ell}^{\mathsf{H}}\mathbf{n}_{k}|^{2}\\}\leq||\mathbf{c}_{k,\ell}||_{2}^{2}\sigma^{2}.$
(79)

Thus, (14) can be then obtained after dividing all terms by
$||\mathbf{c}_{k,\ell}||_{2}^{2}$.

## References

  * [1] W. Jiang, B. Han, M. A. Habibi, and H. D. Schotten, “The road towards 6G: A comprehensive survey,” _IEEE open j. Commun. Soc._ , vol. 2, pp. 334–366, February 2021.
  * [2] A. Magbool, H. Sarieddeen, N. Kouzayha, M.-S. Alouini, and T. Y. Al-Naffouri, “Terahertz-band non-orthogonal multiple access: System- and link-level considerations,” _IEEE Wireless Commun._ , vol. 30, no. 1, pp. 142–149, February 2023.
  * [3] X. Wang _et al._ , “Millimeter wave communication: A comprehensive survey,” _IEEE Commun. Surv. Tutor._ , vol. 20, no. 3, pp. 1616–1653, June 2018.
  * [4] K. Dovelos, S. D. Assimonis, H. Quoc Ngo, B. Bellalta, and M. Matthaiou, “Intelligent reflecting surfaces at terahertz bands: Channel modeling and analysis,” in _2021 IEEE ICC Workshops_ , July 2021, pp. 1–6.
  * [5] E. Basar _et al._ , “Wireless communications through reconfigurable intelligent surfaces,” _IEEE Access_ , vol. 7, pp. 116 753–116 773, August 2019.
  * [6] Z. Zhang _et al._ , “6G wireless networks: Vision, requirements, architecture, and key technologies,” _IEEE Veh. Technol. Mag._ , vol. 14, no. 3, pp. 28–41, July 2019.
  * [7] N. Rajatheva _et al._ , “White paper on broadband connectivity in 6G,” _arXiv:2004.14247_ , April 2020.
  * [8] E. G. Larsson, O. Edfors, F. Tufvesson, and T. L. Marzetta, “Massive MIMO for next generation wireless systems,” _IEEE Commun. Mag._ , vol. 52, no. 2, pp. 186–195, February 2014.
  * [9] T. Huang _et al._ , “A survey on green 6G network: Architecture and technologies,” _IEEE Access_ , vol. 7, pp. 175 758–175 768, December 2019.
  * [10] L. You _et al._ , “Energy efficiency and spectral efficiency tradeoff in RIS-aided multiuser MIMO uplink transmission,” _IEEE Trans. Signal Process._ , vol. 69, pp. 1407–1421, December 2021.
  * [11] M. Zhang, L. Tan, K. Huang, and L. You, “On the trade-off between energy efficiency and spectral efficiency in RIS-aided multi-user MISO downlink,” _Electronics_ , vol. 10, no. 11, May 2021.
  * [12] C. Huang, A. Zappone, G. C. Alexandropoulos, M. Debbah, and C. Yuen, “Reconfigurable intelligent surfaces for energy efficiency in wireless communication,” _IEEE Trans. Wirel. Commun._ , vol. 18, no. 8, pp. 4157–4170, August 2019.
  * [13] Q. Wu _et al._ , “Energy efficient intelligent reflecting surface assisted terahertz communications,” in _2021 IEEE ICC Workshops_ , June 2021, pp. 1–6.
  * [14] Y. Liu, E. Liu, and R. Wang, “Energy efficiency analysis of intelligent reflecting surface system with hardware impairments,” in _2020 IEEE GLOBECOM_ , December 2020, pp. 1–6.
  * [15] A. Mukherjee, V. Kumar, D. W. K. Ng, and L.-N. Tran, “On the energy-efficiency maximization for IRS-assisted MIMOME wiretap channels,” in _2022 IEEE VTC2022-Fall_ , September 2022, pp. 1–6.
  * [16] K. D. Gupta, R. Nigam, D. K. Sharma, and S. K. Dhurandher, “LSTM-based energy-efficient wireless communication with reconfigurable intelligent surfaces,” _IEEE Trans. Green Commun. Netw._ , vol. 6, no. 2, pp. 704–712, June 2022.
  * [17] S. He, Y. Huang, L. Yang, and B. Ottersten, “Coordinated multicell multiuser precoding for maximizing weighted sum energy efficiency,” _IEEE Trans. Signal Process._ , vol. 62, no. 3, pp. 741–751, February 2014.
  * [18] A. Magbool, V. Kumar, and M. F. Flanagan, “On energy efficiency and fairness maximization in RIS-assisted MU-MISO mmwave communications,” in _2023 IEEE ICC 2023_ , June 2023, pp. 5364–5369.
  * [19] C. Lin and G. Y. L. Li, “Terahertz communications: An array-of-subarrays solution,” _IEEE Commun. Mag._ , vol. 54, no. 12, pp. 124–131, December 2016.
  * [20] I. Ahmed _et al._ , “A survey on hybrid beamforming techniques in 5G: Architecture and system model perspectives,” _IIEEE Commun. Surv. Tutor._ , vol. 20, no. 4, pp. 3060–3097, June 2018.
  * [21] G. Zhou, C. Pan, H. Ren, K. Wang, and A. Nallanathan, “A framework of robust transmission design for IRS-aided MISO communications with imperfect cascaded channels,” _IEEE Trans. Signal Process._ , vol. 68, pp. 5092–5106, August 2020.
  * [22] H. Zheng _et al._ , “Robust transmission design for RIS-aided wireless communication with both imperfect CSI and transceiver hardware impairments,” _IEEE Internet Things J._ , vol. 10, no. 5, pp. 4621–4635, November 2023.
  * [23] K. Liu and Z. Zhang, “On the energy-efficiency fairness of reconfigurable intelligent surface-aided cell-free network,” in _2021 IEEE VTC2021-Spring_ , April 2021, pp. 1–6.
  * [24] O. E. Ayach, S. Rajagopal, S. Abu-Surra, Z. Pi, and R. W. Heath, “Spatially sparse precoding in millimeter wave MIMO systems,” _IEEE Trans. Wirel. Commun._ , vol. 13, no. 3, pp. 1499–1513, January 2014.
  * [25] M. R. Akdeniz _et al._ , “Millimeter wave channel modeling and cellular capacity evaluation,” _IEEE J. Sel. Areas Commun._ , vol. 32, no. 6, pp. 1164–1179, June 2014.
  * [26] P. Wang, J. Fang, X. Yuan, Z. Chen, and H. Li, “Intelligent reflecting surface-assisted millimeter wave communications: Joint active and passive precoding design,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 12, pp. 14 960–14 973, October 2020.
  * [27] B. Zheng, C. You, W. Mei, and R. Zhang, “A survey on channel estimation and practical passive beamforming design for intelligent reflecting surface aided wireless communications,” _IEEE Commun. Surv. Tutor._ , vol. 24, no. 2, pp. 1035–1071, February 2022.
  * [28] H. SHI, R. V. Prasad, E. Onur, and I. Niemegeers, “Fairness in wireless networks: Issues, measures and challenges,” _IEEE Commun. Surveys Tuts._ , vol. 16, no. 1, pp. 5–24, May 2014.
  * [29] Q. Shi and M. Hong, “Penalty dual decomposition method for nonsmooth nonconvex optimization—part I: Algorithms and convergence analysis,” _IEEE Trans. Signal Process._ , vol. 68, pp. 4108–4122, June 2020.
  * [30] K. B. Petersen, M. S. Pedersen, _et al._ , “The matrix cookbook,” _Technical University of Denmark_ , vol. 7, no. 15, p. 510, 2008.
  * [31] J. C. Bezdek and R. J. Hathaway, “Convergence of alternating optimization,” _Neural, Parallel & Scientific Computations_, vol. 11, no. 4, pp. 351–368, 2003.
  * [32] A. N. Iusem, “On the convergence properties of the projected gradient method for convex optimization,” _Computational & Applied Mathematics_, vol. 22, pp. 37–52, April 2003.