# Joint Devices and IRSs Association for Terahertz Communications in
Industrial IoT Networks

Muddasir Rahim, Georges Kaddoum, , and Tri Nhu Do M. Rahim, G. Kaddoum, and T.
N. Do are with the Department of Electrical Engineering, the École de
Technologie Supérieure (ÉTS), Université du Québec, Montréal, QC H3C 1K3,
Canada (emails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The Industrial Internet of Things (IIoT) enables industries to build large
interconnected systems utilizing various technologies that require high data
rates. Terahertz (THz) communication is envisioned as a candidate technology
for achieving data rates of several terabits-per-second (Tbps). Despite this,
establishing a reliable communication link at THz frequencies remains a
challenge due to high pathloss and molecular absorption. To overcome these
limitations, this paper proposes using intelligent reconfigurable surfaces
(IRSs) with THz communications to enable future smart factories for the IIoT.
In this paper, we formulate the power allocation and joint IIoT device and IRS
association (JIIA) problem, which is a mixed-integer nonlinear programming
(MINLP) problem. Furthermore, the JIIA problem aims to maximize the sum rate
with imperfect channel state information (CSI). To address this non-
deterministic polynomial-time hard (NP-hard) problem, we decompose the problem
into multiple sub-problems, which we solve iteratively. Specifically, we
propose a Gale-Shapley algorithm-based JIIA solution to obtain stable matching
between uplink and downlink IRSs. We validate the proposed solution by
comparing the Gale-Shapley-based JIIA algorithm with exhaustive search (ES),
greedy search (GS), and random association (RA) with imperfect CSI. The
complexity analysis shows that our algorithm is more efficient than the ES.

###### Index Terms:

Intelligent reconfigurable surfaces (IRSs), Industrial Internet of Things
(IIoT), industrial automation, terahertz (THz).

##  I INTRODUCTION

Following the rapid development of smart manufacturing and industrial
automation, the aim is to increase productivity and efficiency by integrating
the physical world with the internet through the Industrial Internet of Things
(IIoT) [1]. Thousands of industrial devices are now connected to the internet
through the IIoT, enabling them to transmit and gather data that can be used
by business processes and services. However, the increasing demand for
cutting-edge services in innovative industries, like augmented reality
(AR)/virtual reality (VR) maintenance, automatic guided vehicles, and
holographic control display systems, will pose substantial communication
challenges for industrial wireless communications [2]. The large-scale IIoT
networks, being extremely efficient, are intended to provide extremely high-
speed communication of up to terabits-per-second (Tbps) with low latency,
ultra-high reliability, and minimal power consumption [3, 4, 5].

A promising candidate frequency band for wideband communication is terahertz
(THz), which can deliver ultra-high transmission rates, i.e., several Tbps [6,
7]. However, despite the many advantages of the THz band, establishing a
reliable transmission link at THz frequencies remains a challenge in
industrial wireless communications. There are several factors that contribute
to the challenges of establishing a reliable THz transmission link in
industrial wireless communications, including metallic structures, arbitrary
movement of IIoT devices, such as vehicles and robots, factory dimensions, and
electromagnetic interference. Moreover, beam alignment is extremely
challenging due to the narrow beamwidth in THz bands [8, 9, 10].

Intelligent reconfigurable surfaces (IRSs) are gaining increasing attention
from academic as well as industrial communities as a promising technology for
providing smart and reconfigurable industrial wireless communication
environments [11]. The IRS consists of a planar metasurface with multiple
passive reflecting elements. These elements can be digitally controlled in
order to alter the wireless channel between transmitters and receivers by
varying the amplitude and/or phase shift (PS) of the incident signal. This
means IRSs have the potential to reshape the industrial wireless propagation
environment in favor of signal transmission, which is fundamentally different
from existing techniques [12]. Additionally, IRSs have shown potential in
managing dynamic networks, such as dealing with the movement of nodes, changes
in network topology, and network topology reconstruction. IRSs are used to
enhance the coverage area of wireless networks by reflecting and redirecting
the signals to areas that are difficult to reach [13]. This can be especially
useful in urban areas or indoor environments where signals can be obstructed
or weakened by obstacles. IRSs can be used in device-to-device (D2D)
communication to enhance network throughput by creating temporary links
between devices [14]. Furthermore, IRSs can help to mitigate interference
caused by multi-path propagation or co-channel interference. By creating
virtual propagation paths, IRSs can separate signals from different sources,
which can help to improve the signal-to-noise ratio (SNR) and reduce
interference [15].

Moreover, IRSs also enjoy additional practical advantages, such as having a
low profile, lightweight, and conformal geometry, which make them suitable for
wide-scale deployment across industrial wireless networks. Furthermore, an
IRS-assisted network can increase the efficiency of the IIoT network in terms
of data rate, connectivity, and coverage. Considering these facts, we
integrated IRSs into the IIoT networks to facilitate both uplink and downlink
communications. However, optimizing the association time is crucial for
achieving optimal network performance. This is particularly important when the
network topology is subject to frequent changes due to the movement of nodes.
Furthermore, the association time depends on various network conditions,
including node density and mobility. Therefore, it is important to implement
strategies that can reduce the association time to improve the network
performance.

###  I-A Related Works

Recent studies have explored IRSs, potential for improving industrial wireless
communication. Specifically, research works, such as [16, 17, 18, 19, 20, 21],
have studied the energy harvesting, average rate, and reliability analysis of
IRS-assisted industrial networks. For example, S. Dhok et al. integrated IRSs
in an industrial environment to assist the communication between the data
center (DC) and the server machine (SM) [16]. This work aimed to maximize the
system reliability and energy efficiency by optimizing the hybrid power-time
splitting and channel power-based scheduling of destinations. In [17], the
authors investigated the average achievable rate and error probability of IRS-
assisted systems in a factory automation scenario under the finite blocklength
(FBL) regime. Work [17] also derived analytical expressions for both the
average achievable rate and error probability. X. Li et al. considered an IRS-
based wireless-powered non-orthogonal multiple access (NOMA) IoT network to
maximize network sum throughput [18]. The idea behind this work was to use an
alternative optimization (AO) technique to optimize the PS matrices for
wireless energy transfer (WET) and wireless information transfer (WIT).

J. Cheng et al. conducted a study demonstrating the reliability between the
access point (AP) and actuators of multiuser multi-input-single-output (MISO)
in the IRS-assisted IIoT network [19]. The idea of this work was to jointly
optimize the active and passive beamforming, which aims to maximize the number
of actuators with successful decoding for both perfect and imperfect channel
state information (CSI). In [20], the authors demonstrated the potential
features and advantages of IRSs in a smart factory environment to support the
evolution towards Industry $5.0$. Additionally, this work elaborated on the
challenges and potential opportunities for research into future IRS-assisted
wireless factory automation. Finally, M. Schellmann et al. considered an IRS-
aided IIoT communication for a factory environment to maximize the system
capacity [21]. The idea of this work was to derive path-loss models for indoor
factory scenarios with multiple IRSs to assist industrial wireless
communications to maximize the network capacity.

###  I-B Contributions

The aforementioned studies focused on investigating the performance of IRS-
assisted networks by optimizing the beamforming matrix at the IRSs. However,
the joint IIoT device and IRS association (JIIA), a critical factor affecting
industrial wireless communications in THz bands, is largely ignored in the
existing literature. Industrial application demand for high data rates with
good coverage in mobile IIoT device scenarios. Thus, in this work, the IRS is
used to assist the uplink and downlink THz band communication between the IIoT
devices and the AP in an industrial environment. In this context, we formulate
an optimization problem to maximize the sum rate in industrial IoT networks.
This makes it a promising solution for improving the performance of industrial
IoT networks, which are increasingly being used in various industries. The key
contributions of this paper are:

  * •

To the best of the authors’ knowledge, this paper is the first work to study
the power allocation and the JIIA problem in IRS-assisted mobile IIoT networks
operating at THz bands with imperfect CSI. To maximize the sum rate, we
formulate the JIIA problem as a mixed-integer nonlinear programming (MINLP)
problem, which is challenging to solve. Additionally, the IIoT devices are
mobile, which means that the channel condition between IIoT devices and IRSs
may frequently change, requiring a robust algorithm for resource allocation
algorithm including power allocation and IRSs association problem.

  * •

A new framework to tackle the formulated optimization problem. We first
decompose the MINLP problem into multiple sub-problems and then obtain the
optimal PS configuration with a diffuse reflection model. The optimal decoding
and beamforming vectors are achieved using the minimum-mean-square-error
(MMSE) method. Furthermore, we propose an optimal power allocation algorithm
using a water-filling scheme and present the analytical expression for the
optimal power allocation, which facilitates faster algorithm execution.

  * •

We then propose the Gale-Shapley-based robust JIIA algorithm, which obtains a
near-optimal solution to the original problem and significantly reduces the
computational complexity. It is proved that the proposed JIIA algorithm
converges to a stable matching and terminates after a finite number of
iterations in a mobile IIoT scenario.

  * •

Numerical results demonstrate that the proposed JIIA can significantly improve
the sum rate compared to the greedy search (GS) and random association (RA)
approaches. Additionally, the sum rate of our proposed method is comparable to
that of the exhaustive search (ES) method, without considering the association
overhead. This demonstrates the optimality of the proposed JIIA algorithm.
However, when taking into account the association overhead, the sum rate of
the JIIA scheme outperforms the GS, RA, and ES schemes because it minimizes
the negative impact of association overhead.

##  II System Model

Figure 1: Deployment of IRS-assisted THz communications for autonomous robots
in an industrial environment.

We consider IRS-assisted autonomous robots in an IIoT network over THz bands.
As shown in Fig. 1, the access point (AP) is equipped with $K$ antennas for
uplink and downlink communications. Moreover, we assume that the AP operates
in half-duplex (HD) mode and can thus either receive data from the uplink IIoT
device (UD) or transmit data to the downlink IIoT device (DD) at any one
moment. We deployed $I$ uplink IIoT devices and $J$ downlink IIoT devices,
with a single antenna each. In addition, $L$ uplink IRSs and $M$ downlink IRSs
are used to assist the uplink and downlink communications, respectively, where
each IRS is equipped with $N$ reflecting elements. It is assumed that the
direct links from IIoT devices and the AP are unavailable due to unfavorable
propagation conditions111This assumption is reasonable when the direct link is
severely impeded by obstacles or human bodies in the environment [7, 22]..
Therefore, the uplink and downlink communication takes place via IRS deployed
in the network. The area of each IRS element is $A=N_{x}\times N_{y}$, where
$N_{x}$ and $N_{y}$ represent the length of the horizontal and vertical sides,
respectively. In the three-dimensional (3D) Cartesian coordinate system, the
locations of the $i$-th uplink device ($\mathrm{UD}_{i}$) and the $j$-th
downlink device ($\mathrm{DD}_{j}$) at time slot $t$ are denoted by
$\mathbf{c}_{i}(t)=(x_{i}(t),y_{i}(t),z_{i})$ and
$\mathbf{c}_{j}(t)=(x_{j}(t),y_{j}(t),z_{j})$, respectively. Moreover, the
locations of the $l$-th uplink IRS ($\mathrm{UR}_{l}$), the $m$-th downlink
IRS ($\mathrm{DR}_{m}$), and the AP are defined as
$\mathbf{c}_{l}=(x_{l},y_{l},z_{l})$, $\mathbf{c}_{m}=(x_{m},y_{m},z_{m})$,
and
$\mathbf{c}_{\mathrm{AP}}=(x_{\mathrm{AP}},y_{\mathrm{AP}},z_{\mathrm{AP}})$,
respectively. Let $d_{i,l,n}(t)$ denotes the distance from $\mathrm{UD}_{i}$
to the $n$-th element of the $l$-th uplink IRS ($\mathrm{UR}_{l,n}$) at time
slot $t$, and $d_{l,n}$ denotes the distance from $\mathrm{UR}_{l,n}$ to the
AP. For THz communications, the IRS elements are densely deployed, and the
relative distance between the IRS elements is very short. Therefore, we can
ignore the distance between the IRS elements. Thus, the distance $d_{i,l}(t)$
can be computed from the location of $\mathrm{UR}_{l}$ to
$\mathrm{UD}_{i}(t)$, which can be expressed as
$d_{i,l}(t)=\|\mathbf{c}_{i}(t)-\mathbf{c}_{l}\|=\sqrt{(x_{i}(t)-x_{l})^{2}+(y_{i}(t)-y_{l})^{2}+(z_{i}-z_{l})^{2}}.$
Similarly, let $d_{m,n}$ denotes the distance from the AP to the $n$-th
element of the $m$-th downlink IRS ($\mathrm{DR}_{m,n}$), and $d_{j,m,n}(t)$
denotes the distance from $\mathrm{DR}_{m,n}$ to the $j$-th downlink device
$\mathrm{DD}_{j}$ at time slot $t$.

###  II-A Uplink Modeling of the Considered IRS-THz Network

We assume that the direct links between the IIoT devices and the AP are
blocked by an obstacle and thus focus on the propagation model for cascaded
links. Furthermore, we assume that the CSI available at the IIoT devices and
the AP is imperfect222Channel estimation is much more difficult in IRS-
assisted wireless communication than in traditional wireless networks.
Furthermore, it is more challenging at the IRS than the AP because the IRS has
passive reflecting elements that cannot process the pilot signals sent to and
from the IIoT devices / AP [23, 14]. However, channel reciprocity holds for
the IIoT device-to-IRS and IRS-to-AP channels [24]. We therefore use the IRSs’
estimated CSI at the IIoT devices / AP.. Let
$\widehat{\mathbf{g}}_{i,l}\in\mathbb{C}^{N\times 1}$ denote the estimated
channel vector from $\mathrm{UD}_{i}$-to-$\mathrm{UR}_{l}$. In this context,
$[\widehat{\mathbf{g}}_{i,l}]_{n}=\widehat{g}_{i,l,n}$ is the individual
channel from $\mathrm{UD}_{i}(t)$-to-$\mathrm{UR}_{l,n}$, which can be written
as $\widehat{g}_{i,l,n}(t)=\sqrt{\ell_{i,l,n}(t)}e^{-j\omega d_{i,l,n}(t)},$
where $\omega=\frac{2\pi}{\lambda}$ [m-1] is the wavenumber and $\ell_{i,l,n}$
represents the pathloss of the $\mathrm{UD}_{i}$-to-$\mathrm{UR}_{l,n}$ link.
Furthermore, it is important to study the impact of channel estimation errors
(CEEs) to improve optimization design. Then, we define channel with CEEs as
the composition of the estimated channel and its corresponding CEEs, which is
given by

$\displaystyle{g}_{i,l,n}(t)=\widehat{g}_{i,l,n}(t)+\widetilde{g}_{i,l,n}(t),$
(1)

where $\widetilde{g}_{i,l,n}(t)$ denotes the CEEs, which are uncorrelated with
the estimated channel ${g}_{i,l,n}(t)$ and the entries of
$\widetilde{\mathbf{g}}$ are independent and identically distributed (IID)
complex Gaussian with zero mean and variance
$\sigma_{\widetilde{\mathbf{g}}}^{2}$. The correlation between the estimated
and actual channel CSI, which is assumed to be the same for all gains, is
given by

$\displaystyle\rho_{\mathbf{g}}=\frac{\mathbb{E}\\{[\mathbf{g}]_{i}[\widehat{\mathbf{g}}]_{i}^{\ast}\\}}{\sqrt{\mathbb{E}\\{|[\mathbf{g}]_{i}|^{2}\\}\mathbb{E}\\{|[\widehat{\mathbf{g}}]_{i}|^{2}\\}}}=\frac{1}{\sqrt{1+\sigma_{\mathbf{g}}^{2}}},$
(2)

Let $[\widehat{\mathbf{g}}_{l,k}]_{n}=\widehat{g}_{l,n,k}$ be an individual
estimated channel from $\mathrm{UR}_{l,n}$ to antenna $k$ of the AP, where
$\mathbf{g}_{l,k}\in\mathbb{C}^{N\times 1}$ is the channel vector from
$\mathrm{UR}_{l}$ to antenna $k$ of the AP, which can be written as
$\widehat{g}_{l,n,k}=\sqrt{\ell_{l,n,k}}e^{-j\omega d_{l,n,k}},$ where
$\ell_{l,n,k}$ represents the pathloss of the $\mathrm{UR}_{l,n}$-to-antenna
$k$ link. Moreover, with $\widetilde{g}_{l,n,k}$ representing the CEEs, the
channel with CEEs can be written as

$\displaystyle{g}_{l,n,k}=\widehat{g}_{l,n,k}+\widetilde{g}_{l,n,k}.$ (3)

Let $\ell_{i}(t)=\ell_{i,l,n}(t)\ell_{l,n,k}$ denote the overall pathloss of
the cascaded links, which can be expressed as

$\displaystyle\ell_{i}(t)=G_{\mathrm{AP}}G_{{l,n}}(-\mathbf{r}_{i})G_{{l,n}}(\mathbf{r}_{\mathrm{AP}})G_{i}\frac{A^{2}e^{-\kappa_{abs}(f)(d_{i,l}(t)+d_{l})}}{(4\pi
d_{i,l}(t)d_{l})^{2}}.$ (4)

The channel matrix from $\mathrm{UR}_{l}$ to the AP can be denoted as
$\mathbf{G}_{l}=[\mathbf{g}_{l,1},\ldots,\mathbf{g}_{l,k},\ldots,\mathbf{g}_{l,K}]^{\sf
H}\in\mathbb{C}^{K\times N}$. The PS matrix of $\mathrm{UR}_{l}$,
$\mathbf{\Theta}_{l}\in\mathbb{C}^{N\times N}$, is the combination of
amplitude and PS of each element of the $l$-th IRS. The PS matrix of the IRS
can be described as [25]

$\displaystyle\mathbf{\Theta}_{l}=\mathrm{diag}([\kappa_{l,1}e^{j\theta_{l,1}},...,\kappa_{l,n}e^{j\theta_{l,n}},...,\kappa_{l,N}e^{j\theta_{l,N}}]),$
(5)

where $\kappa_{l,n}\in[0,1]$ and $\theta_{l,n}\in[0,2\pi)$ represent the
amplitude and the PS of the $\mathrm{UR}_{l,n}$, respectively. The cascaded
effective channel vector
$\mathbf{h}_{i}^{\mathrm{UR}_{l}}\in\mathbb{C}^{K\times 1}$, from
$\mathrm{UD}_{i}$-to-AP through $\mathrm{UR}_{l}$ can be written as

$\displaystyle\mathbf{h}_{i}^{\mathrm{UR}_{l}}(t)$
$\displaystyle=(\widehat{\mathbf{G}}_{l}+\widetilde{\mathbf{G}}_{l})\mathbf{\Theta}_{l}(\widehat{\mathbf{g}}_{i,l}(t)+\widetilde{\mathbf{g}}_{i,l}(t))=\underbrace{\widehat{\mathbf{G}}_{l}\mathbf{\Theta}_{l}\widehat{\mathbf{g}}_{i,l}(t)}\limits_{\widehat{\mathbf{h}}_{i}^{\mathrm{UR}_{l}}(t)}+\underbrace{\widehat{\mathbf{G}}_{l}\mathbf{\Theta}_{l}\widetilde{\mathbf{g}}_{i,l}(t)+\widetilde{\mathbf{G}}_{l}\mathbf{\Theta}_{l}{\widehat{\mathbf{g}}}_{i,l}(t)}\limits_{\widetilde{\mathbf{h}}_{i}^{\mathrm{UR}_{l}}(t)},$
(6)

where
$\widetilde{\mathbf{h}}_{i}^{\mathrm{UR}_{l}}(t)\sim\mathcal{CN}(\mathbf{0},\sigma_{{\widetilde{\mathbf{g}}}_{i,l}}^{2}\widehat{\mathbf{G}}_{l}\widehat{\mathbf{G}}_{l}^{\sf
H}+\sigma_{\widetilde{\mathbf{G}}_{l}}^{2}\|{\widehat{\mathbf{g}}}_{i,l}\|^{2}\mathbf{I}_{K})\in\mathbb{C}^{K}$.
The channel matrix from the uplink devices to the AP through $\mathrm{UR}_{l}$
can be denoted as
$\mathbf{H}^{\mathrm{UR}_{l}}(t)=[\mathbf{h}_{1}^{\mathrm{UR}_{l}}(t),\ldots,\mathbf{h}_{i}^{\mathrm{UR}_{l}}(t),\ldots,\mathbf{h}_{I}^{\mathrm{UR}_{l}}(t)]\in\mathbb{C}^{K\times
I}$. Considering that $x_{i}$ is a transmitted signal by $\mathrm{UD}_{i}$,
with a transmit power of $p_{i}$. Assuming there are $L$ IRSs, where we use a
single IRS $l$ at a time. The received signals at the AP through
$\mathrm{UR}_{l}$ can be expressed as

$\displaystyle\mathbf{y}_{\mathrm{AP}}^{\mathrm{UR}_{l}}(t)=\textstyle\sum_{i=1}^{I}\sqrt{p_{i}}\widehat{\mathbf{G}}_{l}\mathbf{\Theta}_{l}\widehat{\mathbf{g}}_{i,l}(t)x_{i}(t)+\textstyle\sum_{i=1}^{I}\sqrt{p_{i}}\widehat{\mathbf{G}}_{l}\mathbf{\Theta}_{l}\widetilde{\mathbf{g}}_{i,l}(t)x_{i}(t)+\textstyle\sum_{i=1}^{I}\sqrt{p_{i}}\widetilde{\mathbf{G}}_{l}\mathbf{\Theta}_{l}{\widehat{\mathbf{g}}}_{i,l}(t)x_{i}(t)+\mathbf{n}(t),$
(7)

where $\mathbf{n}\in\mathbb{R}^{K\times 1}$ denotes the AWGN vector at AP, and
$[\mathbf{\Theta}_{l}]_{n}=\kappa_{l,n}e^{j\theta_{l,n}}$ is the $n$-th
element on the diagonal that runs from the top left to the bottom right of the
matrix $\mathbf{\Theta}_{l}$. Moreover, the received signal vector at the AP
can be written in matrix form as

$\displaystyle\mathbf{y}_{\mathrm{AP}}^{\mathrm{UR}_{l}}(t)=\mathbf{H}^{\mathrm{UR}_{l}}(t)\sqrt{\mathbf{P}_{\mathrm{UD}}}\mathbf{x}_{\mathrm{UD}}(t)+\mathbf{n}(t),$
(8)

where $\mathbf{P}_{\mathrm{UD}}\in\mathbb{C}^{I\times I}$ represents the power
matrix of uplink devices.

####  II-A1 Linear detection of the signal at AP

Let $\hat{\mathbf{x}}_{\mathrm{UD}}^{\mathrm{UR}_{l}}(t)\in\mathbb{C}^{I\times
1}$ be the estimated version of $\mathbf{x}_{\mathrm{UD}}(t)$ using a linear
receiver. Considering a linear detection matrix
$\mathbf{V}\in\mathbb{C}^{K\times I}$, which is used to separate the received
signal at the AP into $I$ streams as

$\displaystyle\hat{\mathbf{x}}_{\mathrm{AP}}^{\mathrm{UR}_{l}}(t)$
$\displaystyle=\mathbf{V}^{\sf
H}\mathbf{y}_{\mathrm{AP}}^{\mathrm{UR}_{l}}(t)=\mathbf{V}^{\sf
H}\mathbf{H}^{\mathrm{UR}_{l}}(t)\sqrt{\mathbf{P}_{\mathrm{UD}}}\mathbf{x}_{\mathrm{UD}}(t)+\mathbf{V}^{\sf
H}\mathbf{n}(t).$ (9)

In multi-user scenarios, optimizing channel gains is no longer the sole
objective, as interference must also be considered. While maximizing channel
gains is essential for achieving high signal-to-noise ratios, it is equally
important to minimize interference levels to ensure efficient and reliable
communication among multiple users. This can be achieved through various
interference management techniques, such as beamforming, power control, and
resource allocation. Meanwhile, the $i$-th stream of
$\hat{\mathbf{x}}_{\mathrm{AP}}^{\mathrm{UR}_{l}}(t)$, which is used to detect
$x_{i}$, is given by

$\displaystyle\hat{x}_{\mathrm{AP},i}^{\mathrm{UR}_{l}}(t)$
$\displaystyle={\sqrt{p_{i}}\mathbf{v}_{i}^{\sf
H}\widehat{\mathbf{h}}_{i}^{\mathrm{UR}_{l}}(t)x_{i}(t)}+\sqrt{p_{i}}\mathbf{v}_{i}^{\sf
H}\widetilde{\mathbf{h}}_{i}^{\mathrm{UR}_{l}}(t)x_{i}(t)+{\textstyle\sum_{{i^{\prime}=1,i^{\prime}\neq
i}}^{I}\sqrt{p_{i^{\prime}}}\mathbf{v}_{i}^{\sf
H}{\mathbf{h}}_{i^{\prime}}^{\mathrm{UR}_{l}}x_{i}^{\prime}(t)}+{\mathbf{v}_{i}^{\sf
H}\mathbf{n}(t)},$ (10)

The interference from other devices and noise is treated as effective noise,
and hence, the received signal-to-interference-and-noise ratio (SINR) of
$\mathrm{UD}_{i}$ through $\mathrm{UR}_{l}$ can be expressed as [26]

$\displaystyle\gamma_{i}^{\mathrm{UR}_{l}}(t)=\frac{p_{i}\|\mathbf{v}_{i}^{\sf
H}\widehat{\mathbf{h}}_{i}^{\mathrm{UR}_{l}}(t)|^{2}}{\textstyle\sum_{{i^{\prime}=1,i^{\prime}\neq
i}}^{I}p_{i^{\prime}}|\mathbf{v}_{i}^{\sf
H}\widehat{\mathbf{h}}_{i^{\prime}}^{\mathrm{UR}_{l}}(t)|^{2}+\sum_{i^{\prime}=1}^{I}p_{i^{\prime}}\mathbb{E}\\{|\mathbf{v}_{i}^{\sf
H}\widetilde{\mathbf{h}}_{i^{\prime}}^{\mathrm{UR}_{l}}(t)|^{2}\\}+\sigma_{i}^{2}\|\mathbf{v}_{i}\|^{2}}$
$\displaystyle=\frac{p_{i}\left|\mathbf{v}_{i}^{\sf
H}\widehat{\mathbf{h}}_{i}^{\mathrm{UR}_{l}}(t)\right|^{2}}{\textstyle\sum_{{i^{\prime}=1,i^{\prime}\neq
i}}^{I}p_{i^{\prime}}\left|\mathbf{v}_{i}^{\sf
H}\widehat{\mathbf{h}}_{i^{\prime}}^{\mathrm{UR}_{l}}(t)\right|^{2}+\sum_{i^{\prime}=1}^{I}p_{i^{\prime}}\big{[}\sigma_{\widetilde{\mathbf{g}}_{i^{\prime},l}}^{2}\big{\|}\widehat{\mathbf{G}}_{l}\mathbf{v}_{i}\big{\|}^{2}+\big{(}\sigma_{\widetilde{\mathbf{g}}_{i^{\prime},l}}^{2}\sigma_{\widetilde{\mathbf{G}}_{l}}^{2}+\sigma_{{\widetilde{\mathbf{G}}_{l}}}^{2}\|\widehat{\mathbf{g}}_{i^{\prime},l}\|^{2}\big{)}\|\mathbf{v}_{i}\|^{2}\big{]}+\sigma_{i}^{2}\|\mathbf{v}_{i}\|^{2}}$
(11)

The achievable uplink rate of $\mathrm{UD}_{i}$ through $\mathrm{UR}_{l}$ is
given by

$\displaystyle R_{i}^{\mathrm{UR}_{l}}(t)$
$\displaystyle=\log_{2}(1+{\gamma}_{i}^{\mathrm{UR}_{l}}.$ (12)

The sum rate for the uplink devices through $\mathrm{UR}_{l}$ is the sum of
the individual devices’ rate through $\mathrm{UR}_{l}$, which can be expressed
as

$\displaystyle R_{\tt
sum}^{\mathrm{UR}_{l}}(t)=\textstyle\sum_{i=1}^{I}R_{i}^{\mathrm{UR}_{l}}(t).$
(13)

###  II-B Downlink Modeling of the Considered IRS-THz Network

Consider the downlink channel from the AP with $K$ antennas to a single
antenna downlink device through the $\mathrm{DR}$. Let
$[\widehat{\mathbf{g}}_{k,m}]_{n}=\widehat{g}_{k,m,n}$ be an individual
estimated channel from antenna $k$ of the AP-to-$\mathrm{DR}_{m,n}$, where
$\widehat{\mathbf{g}}_{k,m}\in\mathbb{C}^{N\times 1}$ is the estimated channel
vector from antenna $k$-to-$\mathrm{DR}_{m}$, which can be written as
$\widehat{g}_{k,m,n}=\sqrt{\ell_{k,m,n}}e^{-j\omega d_{m,n}},$ where
$\ell_{k,m,n}$ represents the pathloss over the link from antenna
$k$-to-$\mathrm{DR}_{m,n}$. Furthermore, the channel with CEEs
$\widetilde{g}_{k,m,n}$ can be expressed as

$\displaystyle{g}_{k,m,n}=\widehat{g}_{k,m,n}+\widetilde{g}_{k,m,n}.$ (14)

The channel matrix from the AP-to-$\mathrm{DR}_{m}$ can be denoted as
$\mathbf{G}_{m}=[\mathbf{g}_{m,1},\ldots,\mathbf{g}_{m,k},\ldots,\mathbf{g}_{m,K}]\in\mathbb{C}^{N\times
K}$. The estimated channel from $\mathrm{DR}_{m}$-to-$\mathrm{DD}_{j}$ at time
slot $t$ can be represented as
$\widehat{\mathbf{g}}_{j,m}(t)\in\mathbb{R}^{N\times 1}$. Thus,
$[\widehat{\mathbf{g}}_{j,m}]_{n}=\widehat{g}_{j,m,n}$ is the individual
channel from $\mathrm{DR}_{m,n}$-to-$\mathrm{DD}_{j}$, which can be expressed
as $\widehat{g}_{j,m,n}(t)=\sqrt{\ell_{j,m,n}}e^{-j\omega d_{j,m}},$ where
$\ell_{j,m,n}$ represents the pathloss of the
$\mathrm{UR}_{m,n}$-to-$\mathrm{DD}_{j}$ link. Furthermore, the channel with
CEEs can be written as

$\displaystyle{g}_{j,m,n}(t)=\widehat{g}_{j,m,n}(t)+\widetilde{g}_{j,m,n}(t)$
(15)

where $\widetilde{g}_{j,m,n}(t)$ is an estimation error and is uncorrelated
with the estimated channel ${g}_{j,m,n}(t)$ and the entries of
$\widetilde{\mathbf{g}}$ are IID complex Gaussian with zero mean and variance
$\sigma_{\widetilde{\mathbf{g}}}^{2}$. The correlation between the estimated
and actual channel CSI, which is assumed to be the same for all gains, is
given in (2). Let $\ell_{j}=\ell_{j,m,n}\ell_{k,m,n}$ denote the overall
pathloss of the cascaded downlink, which is given in (4). The cascaded channel
vector $\mathbf{h}_{j}^{\mathrm{DR}_{m}}\in\mathbb{R}^{K\times 1}$ from the AP
to $\mathrm{DD}_{j}$ through $\mathrm{DR}_{m}$ can be expressed as [27]

$\displaystyle\mathbf{h}_{j}^{\mathrm{DR}_{m}}(t)=\underbrace{\widehat{\mathbf{G}}_{m}^{\sf
H}\mathbf{\Theta}_{m}\widehat{\mathbf{g}}_{j,m}(t)}\limits_{\widehat{\mathbf{h}}_{j}^{\mathrm{DR}_{m}}(t)}+\underbrace{\widehat{\mathbf{G}}^{\sf
H}_{m}\boldsymbol{\Theta}_{m}\widetilde{\mathbf{g}}_{j,m}+\widetilde{\mathbf{G}}^{\sf
H}_{m}\boldsymbol{\Theta}_{m}(\widehat{\mathbf{g}}_{j,m}+\widetilde{\mathbf{g}}_{j,m})}\limits_{\widetilde{\mathbf{h}}_{j}^{\mathrm{DR}_{m}}(t)},$
(16)

where
$\widetilde{\mathbf{h}}_{j}^{\mathrm{DR}_{m}}(t)\sim\mathcal{CN}(\mathbf{0},\sigma_{{\widetilde{\mathbf{g}}}_{j,m}}^{2}\widehat{\mathbf{G}}_{m}\widehat{\mathbf{G}}_{m}^{\sf
H}+\sigma_{\widetilde{\mathbf{G}}_{m}}^{2}\|{\widehat{\mathbf{g}}}_{j,m}\|^{2}\mathbf{I}_{K})\in\mathbb{C}^{K}$.
The channel matrix can be written as
$\mathbf{H}^{\mathrm{DR}_{m}}(t)=[\mathbf{h}_{1}^{\mathrm{DR}_{m}}(t),..,\mathbf{h}_{j}^{\mathrm{DR}_{m}}(t),\ldots,\mathbf{h}_{J}^{\mathrm{DR}_{m}}(t)]^{\sf
H}\in\mathbb{C}^{J\times K}$. We consider a precoding matrix
$\mathbf{W}\in\mathbb{C}^{K\times J}$, in which
$\mathbf{W}=[\mathbf{w}_{1},\ldots,\mathbf{w}_{j},\ldots,\mathbf{w}_{J}]$,
where $\mathbf{w}_{j}\in\mathbb{C}^{K\times 1}$ is the precoding vector for
$\mathrm{DD}_{j}$, which can be expressed as
$\mathbf{w}_{j}=\sqrt{p_{j}}\overline{\mathbf{w}}_{j},$ where $p_{j}$ and
$\overline{\mathbf{w}}_{j}$ are the transmit power scaling factor and the
beamforming vector for receiver $j$, respectively. For the given power budget
$\mathbf{P}_{\mathrm{AP}}$, the power constraint can be written as
$\textstyle\sum_{j=1}^{J}\|\mathbf{w}_{j}\|^{2}=\mathsf{Tr}\\{\mathbf{W}\mathbf{W}^{\sf
H}\\}\leq\mathbf{P}_{\mathrm{AP}}.$ Let $\mathbf{x}(t)\in\mathbb{C}^{K\times
1}$ be the transmit signal vector from the AP to $\mathrm{DD}_{J}$, where
$x_{j}(t)$ is the data intended for $\mathrm{DD}_{j}$, then, the transmitted
signal vector from the AP is written as
$\mathbf{x}(t)=\textstyle\sum_{{j=1}}^{J}{\mathbf{w}_{j}}x_{j}(t)=\textstyle\sum_{{j=1}}^{J}\sqrt{p_{j}}\overline{\mathbf{w}}_{j}x_{j}(t)$
[28].

####  II-B1 Linear Precoder at AP

The received signals at $\mathrm{DD}_{j}$ through $\mathrm{DR}_{m}$ at time
slot $t$ can be expressed as

$\displaystyle{y}_{j}^{\mathrm{DR}_{m}}(t)={{\sqrt{p_{j}}(\widehat{\mathbf{h}}_{j}^{\mathrm{DR}_{m}}(t))^{\sf
H}\overline{\mathbf{w}}_{j}x_{{j}}(t)}}+{{\sqrt{p_{j}}(\widetilde{\mathbf{h}}_{j}^{\mathrm{DR}_{m}}(t))^{\sf
H}\overline{\mathbf{w}}_{j}x_{{j}}(t)}}+\\!\\!\\!{\textstyle\sum\limits_{\begin{subarray}{c}j^{\prime}=1,j^{\prime}\neq
j\end{subarray}}^{J}\\!\\!\\!{({\mathbf{h}}_{j}^{\mathrm{DR}_{m}}(t))^{\sf
H}\overline{\mathbf{w}}_{j^{\prime}}x_{j^{\prime}}(t)}\sqrt{p_{j^{\prime}}}}+{n_{j}(t)},$
(17)

where $n_{j}$ is the AWGN with power $\sigma^{2}_{j}$ at $\mathrm{DD}_{j}$.
For a given beamforming matrix333In our model, we consider an IIoT network
with slow mobility, meaning IIoT devices are characterized by limited movement
and relatively stable channel conditions. We therefore use the whole transport
block precoding matrix, which is a suitable assumption.
$\overline{\mathbf{W}}$, the received SINR at $\mathrm{DD}_{j}$ through
$\mathrm{DR}_{m}$ can be expressed as

$\displaystyle\gamma_{j}^{\mathrm{DR}_{m}}(t)=\frac{p_{j}|{(\widehat{\mathbf{h}}_{j}^{\mathrm{DR}_{m}}(t))^{\sf
H}\overline{\mathbf{w}}_{j}}|^{2}}{\sum_{{j^{\prime}=1,j^{\prime}\neq
j}}^{J}p_{j^{\prime}}|{(\widehat{\mathbf{h}}_{j}^{\mathrm{DR}_{m}}(t))^{\sf
H}\overline{\mathbf{w}}_{j^{\prime}}}|^{2}+\sum_{{j^{\prime}=1}}^{J}p_{j^{\prime}}\mathbb{E}\\{|{(\widetilde{\mathbf{h}}_{j}^{\mathrm{DR}_{m}}(t))^{\sf
H}\overline{\mathbf{w}}_{j^{\prime}}}|^{2}\\}+\sigma^{2}_{j}},$
$\displaystyle=\frac{p_{j}\left|{(\widehat{\mathbf{h}}_{j}^{\mathrm{DR}_{m}}(t))^{\sf
H}\overline{\mathbf{w}}_{j}}\right|^{2}}{\sum\limits_{j^{\prime}\neq
j}^{J}p_{j^{\prime}}\left|{(\widehat{\mathbf{h}}_{j}^{\mathrm{DR}_{m}}(t))^{\sf
H}\overline{\mathbf{w}}_{j^{\prime}}}\right|^{2}+\sum\limits_{j^{\prime}=1}^{J}p_{j^{\prime}}\big{[}\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\big{\|}\widehat{\mathbf{G}}_{m}\overline{\mathbf{w}}_{j^{\prime}}\big{\|}^{2}+\big{(}\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\sigma_{\widetilde{\mathbf{G}}_{m}}^{2}+\sigma_{{\widetilde{\mathbf{G}}_{m}}}^{2}\|\widehat{\mathbf{g}}_{j,m}\|^{2}\big{)}\|\overline{\mathbf{w}}_{j^{\prime}}\|^{2}\big{]}+\sigma^{2}_{j}}$
(18)

The achievable downlink rate $R_{j}^{\mathrm{DR}_{m}}$ at $\mathrm{DD}_{j}$
through $\mathrm{DR}_{m}$ can be expressed as

$\displaystyle
R_{j}^{\mathrm{DR}_{m}}(t)=\log_{2}(1+{\gamma}_{j}^{\mathrm{DR}_{m}}(t)).$
(19)

The sum rate of the downlink devices through $\mathrm{UR}_{m}$ is the sum of
the individual devices’ rates, which can be expressed as

$\displaystyle R_{\tt
sum}^{\mathrm{DR}_{m}}(t)=\textstyle\sum_{j=1}^{J}R_{j}^{\mathrm{DR}_{m}}(t).$
(20)

###  II-C Mobility Modeling for The Considered Network

Mobility models are developed to characterize the movement and mobility of
devices as they travel, particularly the changes in their location. These
models are classified into two categories: individual movement models and
group movement models [29]. We consider the random waypoint mobility (RWM) for
moving $\mathrm{UD}_{h}$ and $\mathrm{DD}_{h}$ as a group head, in which
devices move randomly within a specified area. IIoT networks involving freely
moving nodes, such as those on a factory floor or in a warehouse, can be
simulated using this mobility model. Moreover, the remaining $\mathrm{UD}$s
and $\mathrm{DD}$s (group members) follow the respective group head based on
the reference point group mobility (RPGM) [30]. As a result of this model,
nodes travel to a reference point at a fixed speed, stay there for a fixed
amount of time, and then move to the next reference point. Using this mobility
model, one can simulate IIoT networks in which nodes move in coordinated ways,
such as in a manufacturing plant with robots. In the RWM model, the velocity
$\nu_{h}(t)\in[\nu_{min},\nu_{max}]$, where
$h\in\\{\mathrm{UD}_{h},\mathrm{DD}_{h}\\}$ and the direction
$\phi_{h}(t)\in[0,2\pi]$ of the group head at time $t$ are chosen based on a
memoryless random process. In addition, for the given interval of data
collection $T$, we have

$\displaystyle\sin{(\phi_{h}(t))}=\frac{y_{h}(t+T)-y_{h}(t)}{T\nu_{h}(t)},\quad\quad\cos{(\phi_{h}(t))}=\frac{x_{h}(t+T)-x_{h}(t)}{T\nu_{h}(t)}.$
(21)

Figure 2: Device and the AP placement at $t_{1}$ to $t_{4}$ during devices
movement.

The angle deviation factor and the speed deviation factor are denoted as
$\alpha_{a}$ and $\alpha_{s}$, respectively, where $\alpha_{a}\in(-1,1)$ and
$\alpha_{s}\in(-1,1)$. Furthermore, the velocity and the average direction of
group members can be expressed as
$\nu_{{\mathrm{UD}}s/{\mathrm{DD}}s}(t)=|\nu_{h}(t)|+\alpha_{s}\nu_{max}$ and
$\phi_{{\mathrm{UD}}s/{\mathrm{DD}}s}(t)=\phi_{h}(t)+\alpha_{a}\nu_{max}.$ By
incorporating these factors, our model can capture more realistic mobility
patterns of the group members, which is critical for evaluating the
performance of our proposed algorithm. The mobility of devices at $t_{1}$ to
$t_{4}$ is shown in Fig. 2.

##  III Problem Formulation

In IRS-assisted THz networks, there are several unique characteristics and
challenges in terms of network and service management, including limited
coverage, interference, and resource allocation. Furthermore, the network
topology is changing frequently due to the movement of IIoT devices.
Meanwhile, the above-mentioned challenges lead to the need for a robust
resource allocation algorithm including power allocation at the AP and IRS
association algorithm. The achievable rate of the end-to-end transmission
$\mathrm{UR}_{l}\to\mathrm{AP}\to\mathrm{DR}_{m}$ can be expressed as

$\displaystyle R_{l,m}(t)=\big{(}1-\frac{\tau}{T}\big{)}\min\\{R_{\tt
sum}^{\mathrm{UR}_{l}}(t),R_{\tt sum}^{\mathrm{DR}_{m}}(t)\\},$ (22)

where $T$ denotes the coherence interval and $\tau$ denotes the time slots
occupied for JIIA. It is crucial to reduce the association time and employ
proper JIIA techniques to maximize the sum rate. By reducing the association
time, one can allocate more time for data transmission, which can ultimately
improve the overall network sum rate. The reflection optimization at the IRS
and the power allocation optimization at the AP are used to reduce
interference and power consumption, respectively. Additionally, to maximize
the sum rate, proper JIIA techniques must be employed to improve channel gain
and extend the coverage. Let $\mathbf{\Psi}_{L\times M}=[\Psi_{l,m}]$ be the
association matrix, with $\Psi_{l,m}=1$ if $\mathrm{UR}_{l}$ is paired with
$\mathrm{DR}_{m}$ and $\Psi_{l,m}=0$ otherwise. In addition, let
$\mathbf{P}_{1\times J}=[p_{j}]$ be the transmit power vector of the downlink
devices. Then, we can establish the optimization problem as follows.

$\displaystyle(\mathbf{P})\text{ :
}\underset{\mathbf{\Theta_{\varsigma}}\mathbf{V},\mathbf{P},\overline{\mathbf{W}},\mathbf{\Psi}}{\text{maximize}}$
$\displaystyle\textstyle\sum_{l,m}R_{l,m}\Psi_{l,m},$ (23a) subject to
$\displaystyle\textstyle{\theta_{\varsigma,n}}\in[0,2\pi),\forall\varsigma,\forall
n,$ (23b) $\displaystyle\textstyle\kappa_{\varsigma,n}\leq
1,\forall\varsigma,\forall n,$ (23c) $\displaystyle{\|v_{i}\|}=1,\forall{i},$
(23d)
$\displaystyle{\textstyle\sum_{j=1}^{J}p_{j}\|\overline{w}_{j}\|^{2}}\leq\mathbf{P}_{\mathrm{AP}},$
(23e) $\displaystyle\textstyle p_{j}\geq 0,\forall{j},$ (23f)
$\displaystyle\textstyle\sum_{l\in L}\Psi_{l,m}\leq 1,\forall
l\in\\{1,...,L\\},$ (23g) $\displaystyle\textstyle\sum_{m\in M}\Psi_{l,m}\leq
1,\forall m\in\\{1,...,M\\},$ (23h)
$\displaystyle\textstyle\Psi_{l,m}\in\\{0,1\\},\forall l,\forall m,$ (23i)

where $\varsigma\in\\{l,m\\}$. The objective (23a) is the sum rate of the
network with IRS reflection $\mathbf{\Theta}_{\varsigma}$, given decoding
vector $\mathbf{v_{i}}$, beamforming vector $\mathbf{w_{j}}$, downlink devices
transmit power $p_{j}$, and pairing $\Psi_{l,m}$. Constraints (23b) and (23c)
restrict the PS and amplitude of the reflection model, respectively.
Constraints (23e) and (23e) ensure that the sum of the allocated power cannot
exceed the given power budget. Constraint (23g) ensures that each uplink IRS
will be paired to one downlink IRS; constraint (23h) guarantees that each
downlink IRS is allocated to one uplink IRS. Constraints (23g) and (23h)
ensure a one-to-one $\mathrm{UR}\to\mathrm{DR}$ matching. Additionally,
constraint (23i) ensures the matching variables $\Psi_{l,m}$ are binary.

######  Theorem 1.

The optimization problem $(\mathbf{P})$ in (23) is an NP-hard problem.
Furthermore, the problem $(\mathbf{P})$ in (23) is NP-hard, even if we only
consider the IRS association given other variables.

###### Proof:

The proof is given in Appendix A. ∎

Since $\mathbf{P}$ in (23) is NP-hard, solving $\mathbf{P}$ directly becomes
intractable. In the following section, we will reformulate the problem
$\mathbf{P}$ into sub-problems. Then, we develop an optimal solution,
specifically, an exhaustive search-based solution for JIIA. Furthermore, a
low-complexity algorithm based on matching theory will be proposed for JIIA.

###  III-A Problem Decomposition

The formulated problem is an MINLP problem because it includes multiple
integer variables, e.g., $[\mathbf{\Psi}]_{l,m}\in\\{0,1\\}$, and continuous
variables $\theta_{\varsigma,n}\in[0,2\pi]$, and the objective function is
nonlinear [31]. Additionally, the optimization problem in (23) involves
optimizing the IRSs’ reflection, decoding vector, power allocation, precoding
vector, and a matching problem between two disjoint sets (i.e., $\mathrm{UR}$s
and $\mathrm{DR}$s). Thus, our MINLP problem in (23) is NP-hard as it combines
optimization over discrete variables with the challenge of dealing with
nonlinear functions. Therefore its solution usually involves an enormous
search space with exponential time complexity [32].

The goal of (23a) is to maximize the sum rate, while (23b) restricts the PS of
IRS elements, and (23e) ensures that the sum of allocated power cannot exceed
the given power budget. Furthermore, (23g) and (23h) ensure that any
$\mathrm{DR}$ can be assigned to only one $\mathrm{UR}$. This problem is a
variant multiple-knapsack (VMK) problem. The VMK problem is an NP-hard
optimization problem, meaning that it is computationally intractable to find
the optimal solution for large instances of the problem. However, a variety of
heuristics and approximation algorithms have been developed to solve the VMK
problem within a certain tolerance of the optimal solution. In addition, the
VMK problem can be used as a tool for optimization problem decomposition. By
decomposing the optimization problem, it can be more tractable to find a
solution, as the sub-problems can be solved independently and in parallel.
This can also reduce the complexity of the problem, as the sub-problems may be
simpler and easier to solve than the original problem. Hence, we can decompose
this problem into five sub-problems. Firstly, we rewrite the formulated
problem to optimize the PS given the remaining variables as

$\displaystyle(\mathbf{P1})\text{ :
}\underset{{\mathbf{\Theta}_{\varsigma}}}{\text{maximize}}\textstyle(\sum_{l,m}R_{l,m}\Psi_{l,m}|\mathbf{V},\mathbf{P},\overline{\mathbf{W}}),$
$\displaystyle\quad\text{subject
to}\quad\eqref{eq_constraint1a},\eqref{eq_constraint1b}.$ (24)

Then, the optimization of the decoding variable given other variables can be
formulated as

$\displaystyle(\mathbf{P2})\text{ :
}\underset{\mathbf{V}}{\text{maximize}}\textstyle(\sum_{l,m}R_{l,m}\Psi_{l,m}|\mathbf{\Theta}_{\varsigma}^{\star},\mathbf{P},\overline{\mathbf{W}}),$
$\displaystyle\quad\text{subject to}\quad\eqref{eq_constraint2},$ (25)

where $\mathbf{\Theta}_{\varsigma}^{\star}$ is the optimal solution to problem
$(\mathbf{P1})$. In addition, the optimization of downlink devices’ power
allocation given other variables can be expressed as

$\displaystyle(\mathbf{P3})\text{ :
}\underset{\mathbf{P}}{\text{maximize}}\textstyle(\sum_{l,m}R_{l,m}\Psi_{l,m}|\mathbf{\Theta}_{\varsigma}^{\star},\mathbf{V}^{\star},\overline{\mathbf{W}}),$
$\displaystyle\quad\text{subject
to}\quad\eqref{eq_constraint3a},\eqref{eq_constraint3b},$ (26)

where $\mathbf{V}^{\star}$ is the optimal decoding matrix to problem
$(\mathbf{P2})$. Furthermore, we rewrite the optimization problem to optimize
the beamforming vector given other variables as

$\displaystyle(\mathbf{P4})\text{ :
}\underset{\overline{\mathbf{W}}}{\text{maximize}}\textstyle(\sum_{l,m}R_{l,m}\Psi_{l,m}|\mathbf{\Theta}_{\varsigma}^{\star},\mathbf{V}^{\star},\mathbf{P}^{\star}),$
$\displaystyle\quad\text{subject
to}\quad\eqref{eq_constraint3a},\eqref{eq_constraint3b},$ (27)

where ${\mathbf{P}}^{\star}$ is the optimal solution to problem
$(\mathbf{P3})$. Finally, the association of uplink and downlink IRSs given
other variables can be formulated as

$\displaystyle(\mathbf{P5})\text{ :
}\underset{{\mathbf{\Psi}}}{\text{maximize}}\textstyle(\sum_{{l,m}}R_{l,m}\Psi_{l,m}|\mathbf{\Theta}_{\varsigma}^{\star},\mathbf{V}^{\star},\mathbf{P}^{\star},\overline{\mathbf{W}}^{\star}),$
$\displaystyle\quad\text{subject
to}\quad\eqref{eq_constraint4},\eqref{eq_constraint5},\eqref{eq_constraint6},$
(28)

where ${\overline{\mathbf{W}}}^{\star}$ is the optimal solution to problem
$(\mathbf{P4})$.

##  IV Optimal Solutions to Network Design

In this section, we find the optimal solutions for the sub-problems to achieve
a sub-optimal solution for the original problem $(\mathbf{P})$.

###  IV-A Phase Shift Configuration $(\mathbf{P1})$

Reflection methods are classified as specular reflection and diffuse
reflection depending on IRS size and device location. The elements of the IRS
act as perfect mirrors in the specular reflection approach according to
Snell’s reflection law. Meanwhile, the incident wave reflects with the same
beam curvature as the incident beam curvature. However, IRS elements scatter
the incident wave in all directions in the diffuse reflection approach rather
than reflecting with the same beam curvature. For the case of THz
communications, specular reflection is unsuitable due to the absorption and
scattering losses. Therefore, diffuse reflection is utilized in the considered
IIoT networks. Furthermore, we optimize the reflection model to maximize the
achievable data rate with fixed transmit power and decoding matrix/precoding
matrix (24). In order to maximize the received SINR, we need to maximize the
power gain by optimizing the reflection coefficient of IRS elements. By
utilizing (1), (3), and (5), the estimated channel in (6) becomes

$\displaystyle|\widehat{h}_{i}^{\mathrm{UR}_{l}}|^{2}$
$\displaystyle=\textstyle(\sum_{n=1}^{N}\sum_{k=1}^{K}\sqrt{\ell_{i}}e^{-j\omega
d_{l,n,k}}\kappa_{l,n}e^{j\theta_{l,n}}e^{-j\omega
d_{i,l,n}})^{2}=(\sum_{n=1}^{N}\sum_{k=1}^{K}\kappa_{l,n}\sqrt{\ell_{i}}e^{-j\omega(d_{l,n,k}+d_{i,l,n})+j\theta_{l,n}})^{2}.$
(29)

To make the positive contribution of all the reflecting signals, we optimize
the PS by setting $\theta_{l,n}=\omega(d_{l,n,k}+d_{i,l,n})$ and
$\kappa_{l,n}=1$. By doing this, (29) is written as
$|\widehat{h}_{i}^{\mathrm{UR}_{l}}|^{2}=(\textstyle\sum_{n=1}^{N}\sum_{k=1}^{K}\sqrt{\ell_{i}})^{2}.$
We can achieve the optimal solution for the downlink case by following the
same procedure. Furthermore, the optimal PS and amplitude of reflecting
elements at the $\varsigma$-th IRS can be expressed as
$\theta^{\star}_{\varsigma,n}=\omega(d_{i,\varsigma}+d_{\varsigma}),$ and
$\kappa^{\star}_{\varsigma,n}=1,\forall{n},$ respectively. Thus, the optimal
PS configuration at the $\varsigma$-th IRS is determined as

$\displaystyle\mathbf{\Theta}_{\varsigma}^{\star}=\mathrm{diag}[\kappa^{\star}_{\varsigma,1}e^{j\theta^{\star}_{\varsigma,1}},..,\kappa^{\star}_{\varsigma,n}e^{j\theta^{\star}_{\varsigma,n}},..,\kappa^{\star}_{\varsigma,N}e^{j\theta^{\star}_{\varsigma,N}}].$
(30)

###  IV-B Minimum Mean-Square Error (MMSE) Receiver $(\mathbf{P2})$

Here, we find the optimal decoder with fixed transmit power for all uplink
devices and the obtained optimal PS matrix ($\Theta^{\star}$). The MMSE can be
used in a variety of IIoT network applications, including channel estimation,
equalization, and detection. In detection, MMSE can be used to minimize the
mean square error between estimated and actual transmitted signals in the
presence of noise and interference. IIoT networks often involve a large number
of sensors and devices communicating over wireless channels in noisy
environments, making MMSE particularly useful. Therefore, MMSE can improve the
accuracy and reliability of communication systems by reducing the effects of
noise and interference. For the given transmit power, the SINR of all uplink
devices is independently maximized by minimizing the respective MMSE between
the transmitted and estimated signals [33]. As a result, the decoding matrix
can be expressed as

$\displaystyle\mathbf{V}^{\star}$
$\displaystyle=\underset{{\mathbf{V}}}{\mathrm{argmin}}\hskip
1.0pt\mathbb{E}\\{||\mathbf{\hat{x}}_{\mathrm{AP}}^{\mathrm{UR}_{l}}-\mathbf{x}_{\rm
UU}||^{2}\\}=\underset{{\mathbf{V}}}{\mathrm{argmin}}\hskip
1.0pt\mathbb{E}\\{||\mathbf{V}^{\sf
H}\mathbf{y}_{\mathrm{AP}}^{\mathrm{UR}_{l}}-\mathbf{x}_{\rm
UU}||^{2}\\}=\underset{{\mathbf{V}}}{\mathrm{argmin}}\hskip
1.0pt\textstyle\sum_{{\hat{i}=1}}^{I}\mathbb{E}\\{|\mathbf{v}_{i}^{\sf
H}\mathbf{y}_{\mathrm{AP}}^{\mathrm{UR}_{l}}-\mathbf{x}_{i}|^{2}\\}.$ (31)

Moreover, the decoding vector for $\mathrm{UD}_{i}$ can be expressed as

$\displaystyle\mathbf{v}_{i}^{\star}$
$\displaystyle=\underset{{\mathbf{v}_{i}}}{\mathrm{argmin}}\hskip
1.0pt\mathbb{E}\\{|\underbrace{\mathbf{v}_{i}^{\sf
H}\mathbf{y}_{\mathrm{AP}}^{\mathrm{UR}_{l}}}_{\hat{x}_{i}|\mathbf{v}_{i}}-x_{i}|^{2}\\}=\\!\\!(\sigma^{2}_{\widetilde{\mathbf{g}}_{i,l}}\widehat{\mathbf{G}}_{l}\widehat{\mathbf{G}}_{l}^{\sf
H}+\sigma^{2}_{\widetilde{\mathbf{G}}_{l}}\|\widehat{\mathbf{g}}_{i,l}\|^{2}\mathbf{I}_{K})+\sum_{{i^{\prime}=1,i^{\prime}\neq
i}}^{I}p_{i^{\prime}}\mathbf{h}_{i^{\prime},l}\mathbf{h}_{i^{\prime},l}^{\sf
H}\Big{)}^{-1}\widehat{\mathbf{h}}_{i,l},$ (32)

where $\mathbf{I}_{K}$ denotes the $K\times K$ identity matrix. Thus, the SINR
$(\gamma^{\star})$ can be expressed as

$\displaystyle\gamma^{\star}_{\mathrm{UR}_{l},i}=p_{i}\widehat{\mathbf{h}}_{i,l}^{\sf
H}\Big{(}\sigma_{i}^{2}\mathbf{I}_{K}+\sum_{i=1}^{I}p_{i}(\sigma^{2}_{\widetilde{\mathbf{g}}_{i,l}}\widehat{\mathbf{G}}_{l}\widehat{\mathbf{G}}_{l}^{\sf
H}+\sigma^{2}_{\widetilde{\mathbf{G}}_{l}}\|\widehat{\mathbf{g}}_{i,l}\|^{2}\mathbf{I}_{K})+\sum_{{i^{\prime}=1,i^{\prime}\neq
i}}^{I}p_{i^{\prime}}\mathbf{h}_{i^{\prime},l}\mathbf{h}_{i^{\prime},l}^{\sf
H}\Big{)}^{-1}\widehat{\mathbf{h}}_{i,l},$ (33)

###  IV-C Power Allocation Optimization $(\mathbf{P3})$

Power allocation affects IIoT network performance and ultimately impacts the
quality of service (QoS) provided to end users. It is essential to optimize
various network resources, including power allocation, in order to ensure
effective network and service management. Furthermore, power allocation is a
critical resource that must be efficiently managed to maximize network
performance while maintaining QoS standards. In this section, we present a
power allocation optimization algorithm for the AP that uses the Karush-Kuhn-
Tucker (KKT) approach. It is possible to perform this optimization in real-
time, which enables the network to adapt to changing traffic and environmental
conditions in real-time. The optimization problem $(\mathbf{P3})$ to maximize
the sum rate through power allocation can be expressed as

$\displaystyle(\mathbf{P3})\text{ : }$
$\displaystyle\underset{\mathbf{P}}{\text{maximize}}\sum_{j=1}^{J}\log_{2}\Big{(}1+\frac{p_{j}\left|\widehat{h}_{j,j}(t)\right|^{2}}{\textstyle\sum_{{j^{\prime}=1,j^{\prime}\neq
j}}^{J}p_{j^{\prime}}\left|h_{j,j^{\prime}}(t)\right|^{2}+p_{j}\mathcal{X}+\sigma^{2}_{j}}\Big{)},$
(34a) $\displaystyle\text{subject
to}\quad{\textstyle\sum_{{j=1}}^{J}p_{j}\|\overline{w}_{j}\|^{2}}\leq\mathbf{P}_{\mathrm{AP}},\quad\quad\textstyle
p_{j}\geq 0,\quad\quad\forall{j},$ (34b)

where $\widehat{h}_{j,j}(t)$ is the effective channel gain of the desired
device’s signal, ${h}_{j,j^{\prime}}(t)$ is the composition of effective
channel gain and estimation error of the $j^{\prime}$-th device’s signal,
which is an interference for the desired device, and
$\mathcal{X}=\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\|\widehat{\mathbf{G}}_{m}\overline{\mathbf{w}}_{j}\|^{2}+(\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\sigma_{\widetilde{\mathbf{G}}_{m}}^{2}+\sigma_{{\widetilde{\mathbf{G}}_{m}}}^{2}\|\widehat{\mathbf{g}}_{j,m}\|^{2})\|\overline{\mathbf{w}}_{j}\|^{2}$.
The objective of this section is to offer a closed-form solution to the power
allocation problem that supports the offline calculation of power coefficients
using the channel’s CSI. The following theorem presents the closed-form
solution to (34).

######  Theorem 2.

The power allocation matrix is written as $\mathbf{P}_{1\times
J}=\mathsf{diag}\\{p_{1},\ldots,p_{j},\ldots,p_{J}\\}$, and $p_{j}$ is
analytically determined by the water-filling solution as

$\displaystyle
p_{j}=\left[\frac{\frac{1}{\mu\|\overline{\mathbf{w}}_{j}\|^{2}+\upsilon_{j}}-\frac{\iota_{j}(t)+\sigma^{2}_{j}}{\left|\widehat{h}_{j,j}(t)\right|^{2}}}{1+\frac{\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\|\widehat{\mathbf{G}}_{m}\overline{\mathbf{w}}_{j}\|^{2}+(\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\sigma_{\widetilde{\mathbf{G}}_{m}}^{2}+\sigma_{{\widetilde{\mathbf{G}}_{m}}}^{2}\|\widehat{\mathbf{g}}_{j,m}\|^{2})\|\overline{\mathbf{w}}_{j}\|^{2}}{\left|\widehat{h}_{j,j}(t)\right|^{2}}}\right]^{+},$
(35)

where $[x]^{+}=\max\\{x,0\\}$, $\iota_{j}(t)$ is the interference term, and
$\mu$ and $\upsilon_{j}$ are the Lagrangian multipliers. The value of $\mu$ is
selected to satisfy
${\textstyle\sum_{{j=1}}^{J}p_{j}\|\overline{w}_{j}\|^{2}}=\mathbf{P}_{\mathrm{AP}},$
and $\upsilon$ is selected to ensure a non-negative value.

###### Proof:

The proof is given in Appendix B. ∎

We invoke Theorem 2 to obtain the optimal $p_{j}$ and then iterate the water-
filling algorithm over all devices. The iterative water-filling algorithm for
(26) is shown in Algorithm 1.

1 Input:  $J$, power allocation matrix $\mathbf{P}_{1\times J}=[p_{j}]$,
convergence threshold $\epsilon>0$

2 Initialize:  $p_{j}(0)=\mathbf{P}_{AP}/J$, $n=0$

3 for  _$j\in J$_ do

4 while  _( $p_{j}(n)-p_{j}(n-1)>\epsilon$)_ do

5 $\hat{p}_{j}\leftarrow p_{j}(n)$

6 $n=n+1$

7 while  _(until $\hat{p}_{j}$ is not converged)_ do

8 $\iota_{j}=\sum\limits_{j^{\prime}=1,j^{\prime}\neq
j}^{J}\Big{(}p_{j^{\prime}}\left|{(\widehat{\mathbf{h}}_{j}^{\mathrm{DR}_{m}}(t))^{\sf
H}\overline{\mathbf{w}}_{j^{\prime}}}\right|^{2}+\big{[}\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\big{\|}\widehat{\mathbf{G}}_{m}\overline{\mathbf{w}}_{j^{\prime}}\big{\|}^{2}+\big{(}\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\sigma_{\widetilde{\mathbf{G}}_{m}}^{2}+\sigma_{{\widetilde{\mathbf{G}}_{m}}}^{2}\|\widehat{\mathbf{g}}_{j,m}\|^{2}\big{)}\|\overline{\mathbf{w}}_{j^{\prime}}\|^{2}\big{]}\Big{)}$

9 Obtain $\upsilon_{j}$ from equation (13) in [34]

10 Obtain $\mu$ via bisection on
$\mathbf{P}_{\mathrm{AP}}=\sum_{j=1}^{J}\left[{\|\overline{w}_{j}\|^{2}}\Bigg{(}\frac{\frac{1}{\mu\|\overline{w}_{j}\|^{2}+\upsilon_{j}}-\frac{\iota_{j}(t)+\sigma^{2}_{j}}{\left|\widehat{h}_{j,j}(t)\right|^{2}}}{1+\frac{\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\|\widehat{\mathbf{G}}_{m}\overline{\mathbf{w}}_{j}\|^{2}+(\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\sigma_{\widetilde{\mathbf{G}}_{m}}^{2}+\sigma_{{\widetilde{\mathbf{G}}_{m}}}^{2}\|\widehat{\mathbf{g}}_{j,m}\|^{2})\|\overline{\mathbf{w}}_{j}\|^{2}}{\left|\widehat{h}_{j,j}(t)\right|^{2}}}\Bigg{)}\right]^{+}$

11
$p_{j}=\left[\frac{\frac{1}{\mu\|\overline{w}_{j}\|^{2}+\upsilon_{j}}-\frac{\iota_{j}(t)+\sigma^{2}_{j}}{\left|\widehat{h}_{j,j}(t)\right|^{2}}}{1+\frac{\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\|\widehat{\mathbf{G}}_{m}\overline{\mathbf{w}}_{j}\|^{2}+(\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\sigma_{\widetilde{\mathbf{G}}_{m}}^{2}+\sigma_{{\widetilde{\mathbf{G}}_{m}}}^{2}\|\widehat{\mathbf{g}}_{j,m}\|^{2})\|\overline{\mathbf{w}}_{j}\|^{2}}{\left|\widehat{h}_{j,j}(t)\right|^{2}}}\right]^{+}$

12 Normalize to satisfy the power budget

13 $p_{j}(n)\leftarrow\hat{p}_{j}$

14

15

Output:  optimal power allocation matrix $\mathbf{P}^{\star}$

Algorithm 1 Iterative Water-Filling Algorithm for $\mathbf{P3}$ in (26)

###  IV-D Minimum Mean-Square Error (MMSE) based beamforming $(\mathbf{P4})$
for DL transmission

The MMSE beamforming has similar properties as the MMSE receiver; thus, the
SINR at downlink device $j$ is given as

$\displaystyle\gamma_{\mathrm{DR}_{m},j}^{\star}=p_{j}\widehat{\mathbf{h}}_{j,m}^{\sf
H}\left(\sigma_{j}^{2}\mathbf{I}_{K}+\sum_{j=1}^{J}p_{j}(\sigma^{2}_{\widetilde{\mathbf{g}}_{j,m}}\widehat{\mathbf{G}}_{m}\widehat{\mathbf{G}}_{m}^{\sf
H}+\sigma^{2}_{\widetilde{\mathbf{G}}_{m}}\|\widehat{\mathbf{g}}_{j,m}\|^{2}\mathbf{I}_{K})+\sum_{{j^{\prime}=1,j^{\prime}\neq
j}}^{J}p_{j^{\prime}}\mathbf{h}_{j^{\prime},m}\mathbf{h}_{j^{\prime},m}^{\sf
H}\right)^{-1}\widehat{\mathbf{h}}_{j,m}.$ (36)

The MMSE beamforming technique is designed to effectively balance the
mitigation of interference to other users and the Gaussian background noise.
It achieves optimal performance by maximizing the SINR regardless of the SNR
value. Depending on the level of interference, the MMSE beamformer exhibits
characteristics similar to both the zero-forcing beamformer (ZFBF) when
dealing with high inter-user interference, and the maximum ratio transmission
(MRT) beamformer when the interference is low.

###  IV-E Optimal JIIA Strategy using Exhaustive Search (ES) $(\mathbf{P5})$

In this subsection, we focus on the IRS association problem in (23) to find
the optimal association matrix $\mathbf{\Psi}$. In our ES method, the
$\mathrm{UR}\to\mathrm{DR}$ sum rate is utilized to find the optimal pairing
by checking all possible combinations and selecting the optimal one [35, 36].
Thus, our model has $L!$ $\mathrm{UR}$ to $\mathrm{DR}$ combinations. The ES
method for (23) is shown in Algorithm 2. The association matrix and the
achievable rate are set to zero at the beginning of the algorithm.
Furthermore, we check the sum rate of all possible combinations; if the sum
rate of the current combination is greater than the existing combination, we
select the current combination; otherwise, we keep the existing combination.
This process is iterated until we check all combinations, as shown in steps
$4$ to $10$ of Algorithm 2, thus reaching the optimal association matrix
$\mathbf{\Psi}^{\star}$ for JIIA in (28).

1 Input:  $L$, $M$, $R_{l,m}$, permutation of set $L$ as $\mathbf{C}_{\tt
temp}$, association matrix $[\mathbf{\Psi}]_{l,m}$

2 Initialize:  $[\mathbf{\Psi}]_{l,m}=\emptyset$, $\sum R_{l,m}=0$

3 for  _$i\in L!$_ do

4 pairing $\mathrm{UR}_{l}$ with $l$-th element of the $i$-th row in
$\mathbf{C}_{\tt temp}$ as $\mathrm{UR}_{l}\leftarrow\mathbf{C}_{\tt
temp}(i,k)$, $\mathbf{C}_{\tt temp}(i,k)$ gives the $\mathrm{DR}_{m}$

5 for  _$l\in L$_ do

6 $\sum R_{l,m}$ as in equation (22)

7 if ($\sum R_{l,m}\geq\sum R_{l,m}^{\star}$ )

8 $[\mathbf{\Psi}]_{l,m}^{\star}\leftarrow[\mathbf{\Psi}]_{l,m}$

9 else

10 $[\mathbf{\Psi}]_{l,m}^{\star}\leftarrow[\mathbf{\Psi}]_{l,m}^{\star}$

11

Output:  optimal association matrix $[\mathbf{\Psi}]^{\star}_{l,m}$

Algorithm 2 Exhaustive Search-based JIIA Algorithm for $\mathbf{P5}$ in (28)

$(\mathbf{P5})$ in (28).

##  V The Proposed Matching-based JIIA Algorithm

IRS-aided networks require a massive number of reflecting elements $N$,
especially at higher frequencies, where the number of elements can reach
$100\times 100$. The ES method is not an efficient scheme for association
problems. Therefore, to optimize the association problem in (28), we propose a
Gale-Shapley algorithm-based [37] solution. In addition, we assume that the
CSI available at the AP is imperfect. During the signaling phase, priority
matrix configuration and matching are done by the AP based on the CSI and the
JIIA decision is then broadcast to all IIoT devices for the data transmission
phase. The matching algorithm is presented in detail in the following
subsections.

###  V-A Rules for Matching

In this subsection, we introduce some matching rules for the formulated
problem in (28).

######  Definition 1.

(Preference Rule) The $l$-th uplink IRS $\mathrm{UR}_{l}$ prefers the $m$-th
downlink IRS $\mathrm{DR}_{m}$ to the $m^{\star}$-th downlink IRS
$\mathrm{DR}_{m^{\star}}$ is denoted as
$\mathrm{DR}_{m}\succ^{{\mathrm{UR}_{l}}}\mathrm{DR}_{m^{\star}}$.

######  Definition 2.

(Propose Rule) The uplink IRS $l\in L$ proposes to its most preferred downlink
IRS $\mathrm{DR}_{m}$ in its preference list, i.e.,
$\mathrm{DR}_{m}\succ^{{\mathrm{UR}_{l}}}\mathrm{DR}_{m^{\star}}$.

######  Definition 3.

(Reject Rule) The downlink IRS $m\in M$ rejects the proposing uplink IRS if a
better matching uplink IRS exists. Otherwise, the proposed uplink IRS that is
not rejected will be retained as a matching candidate.

######  Definition 4.

(Matching) For the formulated matching problem
$(\mathcal{L},\mathcal{M},[\mathbf{\Upsilon}]_{l,m}^{\mathrm{UR}},[\mathbf{\Upsilon}]_{m,l}^{\mathrm{DR}})$,
where $\mathcal{L}$, $\mathcal{M}$, $[\mathbf{\Upsilon}]_{l,m}^{\mathrm{UR}}$,
and $[\mathbf{\Upsilon}]_{m,l}^{\mathrm{DR}}$ denote the set of uplink IRSs,
the set of downlink IRSs, the priority matrix of uplink IRSs, and the priority
matrix of the downlink IRSs, respectively. The matching $\Psi(m)=j$ indicates
that the downlink IRS $m$ has been matched to uplink IRS $j$, whereas the
matching $\Psi(m)=\emptyset$ implies that the downlink IRS $m$ has not been
allocated to any uplink IRS.

###  V-B Priority Matrix Configuration

The formulated problem, i.e., (28), seeks to maximize the sum rate of the
uplink and downlink through JIIA. The sum rate between $\mathrm{UR}_{l}$ and
$\mathrm{DR}_{m}$ is given in (22). The $\mathrm{UR}$ priorities are obtained
based on the sum rate between each $\mathrm{DR}$ and the $\mathrm{UR}$. After
calculating the sum rate between each $\mathrm{UR}$ and all the downlink IRSs,
the $\mathrm{UR}$ priority matrix is constructed. The priority matrix at
$\mathrm{UR}$s ($[\mathbf{\Upsilon}]_{l,m}^{\mathrm{UR}}$)/$\mathrm{DR}$s
($[\mathbf{\Upsilon}]_{m,l}^{\mathrm{DR}}$) is constructed with the highest
sum rate offered by the $\mathrm{DR}$s/$\mathrm{UR}$s at the top and the
lowest sum rate offered at the bottom, as shown in steps 8 to 19 of Algorithm
3. Additionally, the dimension of the $\mathrm{UR}$ priority matrix
$[\mathbf{\Upsilon}]_{l,m}^{\mathrm{UR}}$ is $L\times M$ and that of the
$\mathrm{DR}$ priority matrix $[\mathbf{\Upsilon}]_{m,l}^{\mathrm{DR}}$ is
$M\times L$. The priority relations for the $\mathrm{UR}$ to the $\mathrm{DR}$
are defined in Definition 1.

###  V-C The Proposed JIIA Algorithm

Input: $L$, $M$, $R_{l,m}$, $[\mathbf{\Upsilon}]_{l,m}^{\mathrm{UR}}$,
$[\mathbf{\Upsilon}]_{m,l}^{\mathrm{DR}}$, set of unmatched $\mathrm{UR}$s
$\Pi$, $[\mathbf{\Psi}]_{l,m}$

1 Initialize $[\mathbf{\Upsilon}]_{l,m}^{\mathrm{UR}}=\emptyset$,
$[\mathbf{\Upsilon}]_{l,m}^{\mathrm{DR}}=\emptyset$,
$[\mathbf{\Psi}]_{l,m}=\emptyset$;

2 Priority Matrix Configuration:

3 for  _$l\in L$_ do

4 $R_{l,m},\quad\forall\hskip 5.69054ptm\in M$ and store in
$[\mathbf{R}]_{l,m}^{\mathrm{UR}}$

5

6$[\mathbf{R}]_{m,l}^{\mathrm{DR}}=[[\mathbf{R}]_{l,m}^{\mathrm{UR}}]^{\sf T}$

7
[$[\mathbf{R}]_{l,m}^{\mathrm{UR}}$,$[\mathbf{\Upsilon}]_{l,m}^{\mathrm{UR}}$]
= sort ($[\mathbf{R}]_{l,m}^{\mathrm{UR}}$,2,descend)

8
[$[\mathbf{R}]_{m,l}^{\mathrm{DR}}$,$[\mathbf{\Upsilon}]_{m,l}^{\mathrm{DR}}$]
= sort ($[\mathbf{R}]_{m,l}^{\mathrm{DR}}$,2,descend)

9 UR - DR Association:

10 while  _(either $\Pi$ $\neq\emptyset$ or $\mathrm{UR}$s not rejected by all
$\mathrm{DR}$s)_ do

11 for  _$\mathrm{UR}_{l^{\prime}}\in\Pi$_ do

12 propose to the highest priority $\mathrm{DR}$ in
$[\mathbf{\Upsilon}]_{l,m}^{\mathrm{UR}}$

13 element of $[\mathbf{\Psi}]_{l,m}=1$

14

15 for  _$\mathrm{DR}_{m}\in M$_ do

16 if ($\mathrm{DR}_{m}\notin[\mathbf{\Psi}]_{l,m}$ )

17
$[\mathbf{\Psi}]_{l,m}\leftarrow[\mathbf{\Psi}]_{l,m}\cup(\mathrm{UR}_{l^{\prime}},\mathrm{DR}_{m})$

18 else if ($R_{m,l^{\prime}}^{\mathrm{DR}}>R_{m,l}^{\mathrm{DR}}$)

19
$[\mathbf{\Psi}]_{l,m}\leftarrow[\mathbf{\Psi}]_{l,m}\cup(\mathrm{UR}_{l^{\prime}},\mathrm{DR}_{m})$

20 else

21
$[\mathbf{\Psi}]_{l,m}\leftarrow[\mathbf{\Psi}]_{l,m}\cup(\mathrm{UR}_{l},\mathrm{DR}_{m})$

22

23

Output:  association matrix $[\mathbf{\overline{\Psi}}]_{l,m}$

Algorithm 3 Proposed Matching-based JIIA Algorithm for $\mathbf{P5}$ in (28)

We deploy a large number of IRSs in the network to assist the communication
and assume that the number of IRSs for uplink communications is the same as
for the downlink (i.e., $L=M$) in the considered network. We explore the JIIA
algorithm as one-to-one matching in the proposed scheme. Each $\mathrm{UR}$
can be allocated to at most one $\mathrm{DR}$ in one-to-one matching.
Algorithm 3 is illustrated as a series of attempts from $\mathrm{UR}$s to
$\mathrm{DR}$s. The $\mathrm{DR}$ can either be associated with the
$\mathrm{UR}$ or remain vacant during this process. Each $\mathrm{UR}$
proposes pairing with the highest priority $\mathrm{DR}$ in their priority
list until the $\mathrm{UR}$ has either been paired or rejected by all
$\mathrm{DR}$s. $\mathrm{UR}$s rejected by the $\mathrm{DR}$ are not allowed
to try again for the same $\mathrm{DR}$. The $\mathrm{DR}$ is immediately
engaged if the $\mathrm{UR}$ proposes a free $\mathrm{DR}$; however, if the
$\mathrm{UR}$ proposes a $\mathrm{DR}$ that is already engaged, the
$\mathrm{DR}$ compares the new and current $\mathrm{UR}$ and gets engaged with
the one with the highest data rate link. Thus, if the $\mathrm{DR}$ favors the
current $\mathrm{UR}$, the new proposal will be rejected. Alternatively, if
the $\mathrm{DR}$ prefers the new proposal, then the $\mathrm{DR}$ will break
the engagement with the current $\mathrm{UR}$ and accept the new
$\mathrm{UR}$. The process is repeated until all $\mathrm{UR}$s are engaged,
or all options are considered, as shown in steps 8 to 19 of Algorithm 3.

######  Remark 1.

The proposed JIIA algorithm terminates in polynomial time, i.e., if there are
$L$ $\mathrm{UR}$ and $M$ $\mathrm{DR}$, then the algorithm is terminated
after at most $L\times M$ iterations.

###### Proof:

In each iteration, each unmatched $\mathrm{UR}$ proposes a $\mathrm{DR}$ that
it has never explored before. Moreover, the $\mathrm{DR}$ accepts or rejects
the new proposal according to the status and preference of the $\mathrm{DR}$,
as shown in steps 8 to 19 of Algorithm 3. As a result, for $L$ $\mathrm{UR}$s
and $M$ $\mathrm{DR}$s, we have $L\times M$ possible proposals occurring in
the proposed algorithm. ∎

######  Theorem 3.

The JIIA matrix $\mathbf{\Psi}^{\star}$ obtained with the ES method and
$\mathbf{\overline{\Psi}}$ obtained with the proposed JIIA algorithm are
similar and converge to the optimal solution for $(\mathbf{P5})$ in (28).

###### Proof:

The proof is given in Appendix C. ∎

###  V-D Computational Complexity and Convergence Analysis

Let’s consider $L=M$ for simplicity and use $L$. Then, the computational
complexity of the ES is $\mathcal{O}(L!)$. Meanwhile, the computational
complexity of the proposed JIIA algorithm is $\mathcal{O}(L^{2})$. To compare
the computational complexity of different methods, we take the natural
logarithm of the linear computational complexity. The logarithm complexity of
the ES is $\mathcal{O}(ln(L!))$, which can be simplified using logarithmic
properties and Stirling’s formula, $\ln(n!)=n\ln(n)-n$. Accordingly, the
logarithm complexity of the ES is $\mathcal{O}(L\ln(L)-L)$. Moreover, the
logarithm complexity of the proposed scheme can be expressed as
$\mathcal{O}(2\ln{L})$, which is much less than that of the ES scheme.

Figure 3: Sum rate versus time slot, for a network area of $40\times 40$
$\mathrm{m}^{2}$, $N=100\times 100$, and $K=64$. Figure 4: Convergence rate
versus number of IRSs, for a network area of $40\times 40$ $\mathrm{m}^{2}$,
$N=100\times 100$, and $K=64$.

Additionally, we consider that the IIoT devices are moving, and the pause time
is short. Thus, the ES may not converge during that short time, while our
proposed scheme can converge. In addition, it can be observed in Fig. 3 that
the proposed scheme converges very fast, while the ES scheme requires more
time slots than the pause time to converge. Furthermore, we examine the
convergence rate with a given pause time for different numbers of IRSs. It is
observed in Fig. 4 that both the ES and proposed schemes converge very fast
for a small number of IRSs. However, when the number of IRSs increases, the
proposed scheme converges in fewer time slots, while the ES requires a huge
time to converge, as shown in Fig. 4.

##  VI Simulation Results and Discussion

In our study, MATLAB is used to implement the proposed algorithm. Each point
is based on the average of $10^{6}$ independent trials. Furthermore, we
present simulation results evaluating the performance of our proposed
algorithm and comparing it to the following schemes as a benchmark.

  * •

Random Association (RA)-based JIIA algorithm: The RA scheme [31] randomly
selects the $\mathrm{UR}-\mathrm{DR}$ association matrix.

  * •

Exhaustive Search (ES)-based JIIA algorithm: In the ES method, the
$\mathrm{UR}-\mathrm{DR}$ sum rate is utilized to find the optimal assignment
by checking all possible combinations and selecting the optimal one, as in
Section IV-E.

  * •

Greedy Search (GS)–based JIIA algorithm: In the GS method [38], each
$\mathrm{UR}$ selects the $\mathrm{DR}$ with the highest data rate.
$\mathrm{DR}$s receiving multiple proposals are randomly allocated to one of
the $\mathrm{UR}$s.

###  VI-A Simulation Setup

In the simulated scenario, explained in Section II, $10$ uplink and downlink
devices are deployed within the network area. Additionally, we deployed IRSs
with $100\times 100$ elements at each IRS. The length and width of each IRS
element are $N_{x}=N_{y}=0.4\lambda$ [39], where $\lambda$ is the wavelength.
Table I details the simulation parameters based on existing works. In our
study, we used reference point group mobility mode (RPGM) with synchronized
velocity and pause time of all devices; however, the directions are different,
as shown in Fig 5. The snapshot of a single time slot is shown in Fig 5 (a),
i.e., the state of the devices and their positions at a specific moment in
time. Furthermore, Fig 5 (b) depicts the mobility of devices in different
directions.

TABLE I: Simulation Parameters Parameters | Values |  Parameters | Values
---|---|---|---
Carrier frequency, $f_{c}$ [GHz] | $300$ |  Number of AP antennas, $K$ | $64$
Number of $\mathrm{DD}$, $J$ | $10$ |  Absorption coefficient, $\kappa_{abs}(f)$ ${[}m^{-1}{]}$ | $0.0033$
Number of $\mathrm{UD}$, $I$ | $10$ |  Noise power density, $N_{0}$ [dBm/Hz] | $-174$
Transmit power, $P$ [dBm] | $23$ |  Network area [$\mathrm{m}^{2}$] | $40\times 40$
Coherence interval | $200$ |  Number of IRS elements, $N$ | $100\times 100$
Channel bandwidth, $B$ [GHz] | $10$ |  Noise figure, NF [dB] | $10$
Spacing between IRS elements | $\lambda/2$ |  Angle deviation factor, $\alpha_{a}$ | $(-1,1)$
Spacing between AP antennas | $\lambda/2$ |  Speed deviation factor, $\alpha_{s}$ | $(-1,1)$
Side length of IRS elements | $0.4\lambda$ |  $\mathrm{UD}_{i}$ coordinate $[m]$ |  $x_{i}\sim\mathcal{U}(0\,\,20)$, $y_{i}\sim\mathcal{U}(0\,\,40)$
Amplitude of IRS element coef., $\kappa$ | $1$ |  $\mathrm{DD}_{j}$ coordinate $[m]$ |  $x_{j}\sim\mathcal{U}(20\,\,40)$, $y_{j}\sim\mathcal{U}(0\,\,40)$
PS of IRS elements, $\theta$ | $[0,2\pi]$ |  $\mathrm{UR}_{l}$ coordinate $[m]$ |  $x_{l}\sim\mathcal{U}(5\,\,20)$, $y_{l}\sim\mathcal{U}(0\,\,40)$
Avg. height of IRSs’ location $[m]$ | $10$ |  $\mathrm{DR}_{m}$ coordinate $[m]$ |  $x_{m}\sim\mathcal{U}(20\,\,35)$, $y_{m}\sim\mathcal{U}(0\,\,40)$
Avg. height of IIoT devices $[m]$ | $1$ |  $\mathrm{AP}$ coordinate $[m]$ | $[20,20,10]$
Figure 5: Simulation scenario for a network area of $40\times 40$
$\mathrm{m}^{2}$ and $I=J=10$; (a) single snapshot of device deployment, (b)
device deployment with mobility.

###  VI-B Performance Comparison without Association Overhead

This section evaluates the performance of the proposed algorithm based on the
network sum rate without considering the association overhead (i.e.,
$\tau=0$). The results demonstrate how transmission power, number of antennas
at AP, number of IRS elements, and size of the network impact the sum rate.

####  VI-B1 Impact of the transmission power

To demonstrate how transmission power affects performance, we use $I=J=10$,
$K=64$, $N=100\times 100$, and a network area of $40\times 40$
$\mathrm{m}^{2}$. Fig. 6 shows the sum rate under different transmission
powers. For all cases, the sum rate increases with the transmission power.
Additionally, the sum rate of the proposed scheme matches that of the ES
scheme. Furthermore, the proposed algorithm outperforms the GS and RA
algorithms.

Figure 6: Sum rate versus transmission power, for a network area of $40\times
40$ $\mathrm{m}^{2}$, $N=100\times 100$, $I=J=10$, and $K=64$.

####  VI-B2 Impact of the number of antennas at the AP

Figure 7: Sum rate versus number of antennas at the AP, for a network area of
$40\times 40$ $\mathrm{m}^{2}$, $N=100\times 100$, and $P=23$ dBm.

Comparing the performance of conventional large-scale antenna arrays deployed
at the AP is imperative. The sum rate increases as the number of antennas at
the AP increases, as shown in Fig. 7. This is because when $K$ increases, the
total signal gain increases as well, which increases the sum rate.
Furthermore, it is observed from Fig. 7 that a significant sum rate
performance improvement can be achieved by the proposed scheme over the GS and
RA schemes.

####  VI-B3 Impact of the number of IRS elements

In Fig. 8, the sum rate is plotted against the number of IRS elements.
According to the results, the sum rate of all schemes increases with $N$. This
is because the received signal strength increases as $N$ increases, which
leads to an improvement in the sum rate. Moreover, the results show that the
sum rate of the proposed scheme is similar to that of the ES scheme.
Additionally, the proposed scheme achieves a higher sum rate than the GS and
RA schemes, due to proper association.

Figure 8: Sum rate versus number of IRS elements, for a network area of
$40\times 40$ $\mathrm{m}^{2}$, $K=64$, $I=J=10$, $K=64$, and $P=23$ dBm.

####  VI-B4 Impact of the carrier frequency

In Fig. 9, the sum rate is plotted against the number of IRS elements for
three different frequencies. The sum rate of all the schemes increases with
the number of IRS elements in all three cases. Furthermore, our proposed
method achieves a higher sum rate than the GS and RA schemes in all three
frequencies. The main purpose of Fig. 9 is to compare the sum rates at
different frequencies. We therefore measured the sum rate in bits per second
(bps) instead of bps/Hz to show the improvement for the THz band. For a fixed
number of IRS elements, i.e., $100\times 100$, the sum rates achieved at $6$
GHz, $28$ GHz, and $300$ GHz are $0.55$ Gbps, $2.4$ Gbps, and $4$ Gbps,
respectively, as shown in Fig. 9. Thus, the sum rate at the proposed frequency
(i.e., $300$ GHz) is approximately $40\%$ higher than at $28$ GHz and $87\%$
higher than at $6$ GHz.

(a) $f_{c}$= $6$ GHz

(b) $f_{c}$= $28$ GHz (c) $f_{c}$= $300$ GHz

Figure 9: Sum rate versus number of IRS elements, for a network area of
$40\times 40$ $\mathrm{m}^{2}$, $K=64$, and $P=23$ dBm.

####  VI-B5 Impact of the Network Area

Fig. 10 shows that the sum rate decreases with increases in the network size
for all schemes due to increased distance between devices. In particular, the
pathloss is distance-dependent, so as distance increases, the sum rate
decreases. In addition, the sum rate achieved by the proposed scheme is
approximately identical to that of the ES scheme and superior to the sum rates
of the GS and RA schemes, as shown in Fig. 10.

Figure 10: Sum rate versus network area, for $N=100\times 100$, $K=64$,
$I=J=10$, and $P=23$ dBm.

####  VI-B6 Impact of the Mobility

Figure 11: Sum rate versus location index for different mobility scenarios,
for a network area of $40\times 40$ $\mathrm{m}^{2}$, $K=64$, and $P=23$ dBm.

In Fig. 11, the performance of the proposed JIIA and the ES algorithms is
compared in different mobility scenarios. The JIIA algorithm achieves a
similar sum rate in both slow and fast mobility scenarios, indicating that it
is robust to changes in mobility. In contrast, the ES algorithm performs
comparably to the JIIA algorithm only in scenarios with slow mobility.
However, in scenarios with fast mobility, the ES algorithm achieves a lower
sum rate than the proposed JIIA algorithm. This is because the ES algorithm
requires a longer period to search for the optimal association, while the
pause time is short in fast mobility scenarios.

####  VI-B7 Impact of the CSI

In this section, we demonstrate the impact imperfect CSI has and compare the
results obtained with perfect CSI and imperfect CSI. Fig. 12 shows that the
sum rate is lower with imperfect CSI for all the schemes. This is due to
channel uncertainty. Furthermore, it can be observed that the proposed scheme
performs similarly to the ES method in both the perfect and imperfect CSI
cases.

Figure 12: Sum rate versus transmission power, for a network area of $40\times
40$ $\mathrm{m}^{2}$, $N=100\times 100$, $I=J=10$, and $K=64$.

In Fig. 13, the sum rate is plotted against the CSI condition parameter. The
results indicate the sum rate decreases with the CSI condition parameter for
all the schemes. This is due to changes in the channel condition. Moreover,
the results show that the sum rate of the proposed scheme is approximately
similar to that of the ES method. Furthermore, the proposed scheme achieves a
higher sum rate than the GS and RA schemes.

Figure 13: Sum rate versus CSI condition parameter, for a network area of
$40\times 40$ $\mathrm{m}^{2}$, $N=100\times 100$, $I=J=10$, and $K=64$.

###  VI-C Performance Comparison with Association Overhead

The objective of this section is to assess the effectiveness of the proposed
algorithm in terms of network sum rate while taking into account association
overhead. To estimate association overhead, we count the number of time slots
required for JIIA and then divide it by the coherence interval. Furthermore,
by subtracting this overhead duration from 1 and multiplying that number by
the achievable rate, we obtain the achievable rate (22) with association
overhead.

Figure 14: Sum rate versus transmission power considering association
overhead, for a network area of $40\times 40$ $\mathrm{m}^{2}$, $N=100\times
100$, $I=J=10$, and $K=64$.

The results presented in Fig. 14 demonstrate the relationship between
transmission power and sum rate, as well as the comparative performance of
different algorithms. Across all considered cases, we observe that increasing
the transmission power results in a corresponding increase in the sum rate.
Moreover, the proposed algorithm achieves superior performance compared to the
ES, GS, and RA schemes. This performance advantage is attributed to the
proposed algorithm’s ability to minimize the negative impact of the
association overhead, which is a significant drawback of the ES scheme.

##  VII Conclusions

This paper investigated IRS-assisted THz communications to enable future smart
industries for the IIoT, where several IRSs are deployed to assist with
communication. First, we formulated the channel model with perfect CSI and
imperfect CSI. Then, we formulated an optimization problem to maximize the
network sum rates. It was NP-hard, so we decomposed it to solve it. We solved
the PS configuration with a diffuse reflection model. Then, the decoding and
beamforming vectors are achieved using the MMSE method. In this work, we aimed
to solve the power allocation problem $(\mathbf{P3})$ and the JIIA problem
$(\mathbf{P5})$. We achieved the analytical expression for the optimal power
allocation for $(\mathbf{P3})$. Furthermore, we solved the JIIA problem
$(\mathbf{P5})$ with our proposed matching-based JIIA algorithm. Simulation
results validated the potential of the proposed IRS-assisted THz communication
scheme and demonstrated its superiority over GS and RA schemes. We also
analyzed the performance while considering the association overhead, and the
results showed that the proposed scheme outperforms the ES scheme. This is
because the ES scheme incurs a high association overhead, whereas our proposed
scheme does not. Additionally, the computational complexity analysis confirmed
that the proposed JIIA algorithm is more computationally efficient than the
ES.

##  Appendix A Proof of Theorem 1

To demonstrate that problem (9) is NP-hard using computational complexity
theory, we can follow these three steps: (1) Select a well-known decision
problem $Q$ that is already proven to be NP-complete. (2) Create a polynomial
time conversion from any instance of $Q$ to an instance of problem (23). (3)
Prove that both instances have the same objective value after the
transformation. We first assume the case of power allocation and JIIA with the
given remaining variables. The sum-rate maximization problem in (23) can be
expressed as

$\displaystyle\underset{\mathbf{P},\mathbf{\Psi}}{\text{maximize}}$
$\displaystyle\textstyle\sum_{l,m}R_{l,m}\Psi_{l,m},$ (37a) subject to
$\displaystyle{\textstyle\sum_{j=1}^{J}p_{j}\|\overline{w}_{j}\|^{2}}\leq\mathbf{P}_{\mathrm{AP}},\quad\textstyle
p_{j}\geq 0,\forall{j},$ (37b) $\displaystyle\textstyle\sum_{l\in
L}\Psi_{l,m}\leq 1,\forall{l\in L},\quad\textstyle\sum_{m\in M}\Psi_{l,m}\leq
1,\forall{m\in M},$ (37c)
$\displaystyle\textstyle\Psi_{l,m}\in\\{0,1\\},\forall{l,m},$ (37d)

which has been proved to be NP-hard in [31]. Furthermore, we prove that the
IRS association is NP-hard for the given power allocation matrix. The sum-rate
maximization problem in (23) can be expressed as

$\displaystyle\underset{\mathbf{\Psi}}{\text{maximize}}$
$\displaystyle\textstyle\sum_{l,m}R_{l,m}\Psi_{l,m},\quad\text{subject
to}\quad\eqref{eq_constraint55},\eqref{eq_constraint66},$ (38)

which is similar to the link scheduling problem, a known NP-hard problem [40].
This completes the proof of Theorem 1.

##  Appendix B Proof of Theorem 2

The Lagrangian function associated with problem $(\mathbf{P3})$ in (34) can be
expressed as

$\displaystyle
L(p_{j},\mu)=\sum\limits_{\begin{subarray}{c}j=1\end{subarray}}^{J}\log_{2}\bigg{(}1+\frac{p_{j}\left|\widehat{h}_{j,j}(t)\right|^{2}}{\textstyle\sum_{{j^{\prime}=1,j^{\prime}\neq
j}}^{J}p_{j^{\prime}}\left|h_{j,j^{\prime}}(t)\right|^{2}+p_{j}\mathcal{X}+\sigma^{2}_{j}}\bigg{)}-\mu\big{(}\textstyle\sum_{{j=1}}^{J}p_{j}\|\overline{w}_{j}\|^{2}-P_{\mathrm{AP}}\big{)}-\upsilon_{j}p_{j},$
(39)

where
$\mathcal{X}=\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\|\widehat{\mathbf{G}}_{m}\overline{\mathbf{w}}_{j}\|^{2}+(\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\sigma_{\widetilde{\mathbf{G}}_{m}}^{2}+\sigma_{{\widetilde{\mathbf{G}}_{m}}}^{2}\|\widehat{\mathbf{g}}_{j,m}\|^{2})\|\overline{\mathbf{w}}_{j}\|^{2}$.
By setting the derivative of (39) with respect to $p_{j}$ to zero, we can
obtain the KKT system of the optimization problem as

$\displaystyle\frac{\partial L(p_{j},\mu)}{\partial p_{j}}=0$ $\displaystyle
p_{j}\Big{(}1+\frac{\mathcal{X}}{\left|\widehat{h}_{j,j}(t)\right|^{2}}\Big{)}+\frac{\sum_{{j^{\prime}=1,j^{\prime}\neq
j}}^{J}p_{j^{\prime}}\left|h_{j,j^{\prime}}(t)\right|^{2}+\sigma^{2}_{j}}{\left|\widehat{h}_{j,j}(t)\right|^{2}}=\frac{1}{\mu\|\overline{w}_{j}\|^{2}+\upsilon_{j}},$
(40)

where $\mu$ and $\upsilon_{j}$ are the Lagrangian multipliers. The multiplier
$\upsilon_{j}$ stands for the effect of interference caused by the $j$-th
device on other devices and is called the taxation term [34]. The interference
term can be represented as

$\displaystyle\iota_{j}(t)$
$\displaystyle=\textstyle\sum_{{j^{\prime}=1,j^{\prime}\neq
j}}^{J}p_{j^{\prime}}\left|h_{j,j^{\prime}}(t)\right|^{2},$
$\displaystyle=\sum\limits_{j^{\prime}=1,j^{\prime}\neq
j}^{J}\Big{(}p_{j^{\prime}}\left|{(\widehat{\mathbf{h}}_{j}^{\mathrm{DR}_{m}}(t))^{\sf
H}\overline{\mathbf{w}}_{j^{\prime}}}\right|^{2}+\big{[}\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\big{\|}\widehat{\mathbf{G}}_{m}\overline{\mathbf{w}}_{j^{\prime}}\big{\|}^{2}+\big{(}\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\sigma_{\widetilde{\mathbf{G}}_{m}}^{2}+\sigma_{{\widetilde{\mathbf{G}}_{m}}}^{2}\|\widehat{\mathbf{g}}_{j,m}\|^{2}\big{)}\|\overline{\mathbf{w}}_{j^{\prime}}\|^{2}\big{]}\Big{)}.$
(41)

To find the power of the $j$-th device, we rearrange (B) and utilize (B) as

$\displaystyle{p_{j}\Big{(}1+\frac{\mathcal{X}}{\left|\widehat{h}_{j,j}(t)\right|^{2}}\Big{)}+\frac{\textstyle\sum_{{j^{\prime}=1,j^{\prime}\neq
j}}^{J}p_{j^{\prime}}|h_{j,j^{\prime}}(t)|^{2}}{|\widehat{h}_{j,j}(t)|^{2}}+\frac{\sigma^{2}_{j}}{|\widehat{h}_{j,j}(t)|^{2}}}=\frac{1}{\mu\|\overline{w}_{j}\|^{2}+\upsilon_{j}}$
$\displaystyle
p_{j}\Big{(}1+\frac{\mathcal{X}}{|\widehat{h}_{j,j}(t)|^{2}}\Big{)}=\frac{1}{\mu\|\overline{w}_{j}\|^{2}+\upsilon_{j}}-\frac{\iota_{j}(t)+\sigma^{2}_{j}}{|\widehat{h}_{j,j}(t)|^{2}}$
$\displaystyle
p_{j}=\left[\frac{\frac{1}{\mu\|\overline{w}_{j}\|^{2}+\upsilon_{j}}-\frac{\iota_{j}(t)+\sigma^{2}_{j}}{\left|\widehat{h}_{j,j}(t)\right|^{2}}}{1+\frac{\mathcal{X}}{\left|\widehat{h}_{j,j}(t)\right|^{2}}}\right]^{+}=\left[\frac{\frac{1}{\mu\|\overline{w}_{j}\|^{2}+\upsilon_{j}}-\frac{\iota_{j}(t)+\sigma^{2}_{j}}{\left|\widehat{h}_{j,j}(t)\right|^{2}}}{1+\frac{\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\|\widehat{\mathbf{G}}_{m}\overline{\mathbf{w}}_{j}\|^{2}+(\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\sigma_{\widetilde{\mathbf{G}}_{m}}^{2}+\sigma_{{\widetilde{\mathbf{G}}_{m}}}^{2}\|\widehat{\mathbf{g}}_{j,m}\|^{2})\|\overline{\mathbf{w}}_{j}\|^{2}}{\left|\widehat{h}_{j,j}(t)\right|^{2}}}\right]^{+},$
(42)

where $[x]^{+}=\max\\{x,0\\}$. Given $\iota_{j}$ and $\upsilon_{j}$, the
water-filling algorithm can find the optimal $p_{j}$. Substituting (42) into
constraint (34b) yields

$\displaystyle\mathbf{P}_{\mathrm{AP}}=\sum_{j=1}^{J}\left[{\|\overline{w}_{j}\|^{2}}\Bigg{(}\frac{\frac{1}{\mu\|\overline{w}_{j}\|^{2}+\upsilon_{j}}-\frac{\iota_{j}(t)+\sigma^{2}_{j}}{\left|\widehat{h}_{j,j}(t)\right|^{2}}}{1+\frac{\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\|\widehat{\mathbf{G}}_{m}\overline{\mathbf{w}}_{j}\|^{2}+(\sigma_{\widetilde{\mathbf{g}}_{j,m}}^{2}\sigma_{\widetilde{\mathbf{G}}_{m}}^{2}+\sigma_{{\widetilde{\mathbf{G}}_{m}}}^{2}\|\widehat{\mathbf{g}}_{j,m}\|^{2})\|\overline{\mathbf{w}}_{j}\|^{2}}{\left|\widehat{h}_{j,j}(t)\right|^{2}}}\Bigg{)}\right]^{+}.$
(43)

This completes the proof of Theorem 2.

##  Appendix C Proof of Theorem 3

The MINLP problem formulated is NP-hard since it combines optimizing over
discrete variables with nonlinear functions. In this context, the ES scheme is
conventionally used to find the optimal solution [36]. The optimality of the
ES method can be proven by the following conditions.

  * •

The ES method’s characteristics dictate that, at every iteration, the sum rate
for the current association matrix is calculated as in step 7 of Algorithm 2.

  * •

Then, in steps 9 to 12 of Algorithm 2, we check if the sum rate of the current
association matrix is greater than that of the previous one and select the
matrix with the higher sum rate.

This way, we obtain the optimal association matrix $\mathbf{\Psi}^{\star}$,
resulting in the highest sum rate.

Now, we prove that the JIIA matrix achieved with the proposed algorithm is
stable and optimal. First, we note that the algorithm must terminate in
polynomial time, as indicated in Remark 1. When it terminates, we have a
match. Let’s prove it is stable and optimal. Suppose the proposed algorithm
matches $\mathrm{UR}_{l}$ and $\mathrm{DR}_{m}$. However, $\mathrm{UR}_{l}$
prefers $\mathrm{DR}_{m^{\prime}}$ over $\mathrm{DR}_{m}$. When the algorithm
is run, $\mathrm{UR}_{l}$ offers to $\mathrm{DR}$s in order of preference, so
if $\mathrm{UR}_{l}$ ends up with $\mathrm{DR}_{m}$, it must have proposed to
$\mathrm{DR}_{m^{\prime}}$ at some point and been rejected. Now, when
$\mathrm{DR}_{m^{\prime}}$ rejected $\mathrm{UR}_{l}$, it must have had
another offer in hand that it preferred over that of $\mathrm{UR}_{l}$. As the
algorithm proceeded, $\mathrm{DR}_{m^{\prime}}$ must have received a better
offer, so it ended up with a $\mathrm{UR}$ that is preferred over
$\mathrm{UR}_{l}$. Therefore, $\mathrm{DR}_{m^{\prime}}$ will not form a
blocking pair with $\mathrm{UR}_{l}$. It follows that, at the end of the
algorithm, there will be no better uplink IRS $\mathrm{UR}_{l}$ that can form
a blocking pair with downlink IRS $\mathrm{DR}_{m}$ than its current partner.
Hence the matching is stable and optimal. This completes the proof of Theorem
3.

## References

  * [1] E. Sisinni, A. Saifullah, S. Han, U. Jennehag, and M. Gidlund, “Industrial internet of things: Challenges, opportunities, and directions,” _IEEE Trans. Ind. Inf._ , vol. 14, no. 11, pp. 4724–4734, 2018.
  * [2] K. B. Letaief, W. Chen, Y. Shi, J. Zhang, and Y.-J. A. Zhang, “The roadmap to 6G: AI empowered wireless networks,” _IEEE Commun. Mag._ , vol. 57, no. 8, pp. 84–90, 2019.
  * [3] S. Zeb, A. Mahmood, S. A. Hassan, M. J. Piran, M. Gidlund, and M. Guizani, “Industrial digital twins at the nexus of nextG wireless networks and computational intelligence: A survey,” _Journal of Network and Computer Applications_ , p. 103309, 2022.
  * [4] N. Rubab, S. Zeb, A. Mahmood, S. A. Hassan, M. I. Ashraf, and M. Gidlund, “Interference mitigation in RIS-assisted 6G systems for indoor industrial IoT networks,” in _2022 IEEE 12th Sensor Array and Multichannel Signal Processing Workshop (SAM)_. IEEE, 2022, pp. 211–215.
  * [5] W. Mao, Z. Zhao, Z. Chang, G. Min, and W. Gao, “Energy-efficient industrial internet of things: overview and open issues,” _IEEE Trans. Ind. Inf._ , vol. 17, no. 11, pp. 7225–7237, 2021.
  * [6] Z. Chen, X. Ma, B. Zhang, Y. Zhang, Z. Niu, N. Kuang, W. Chen, L. Li, and S. Li, “A survey on terahertz communications,” _China Communications_ , vol. 16, no. 2, pp. 1–35, 2019.
  * [7] Z. Wan, Z. Gao, F. Gao, M. Di Renzo, and M.-S. Alouini, “Terahertz massive MIMO with holographic reconfigurable intelligent surfaces,” _IEEE Trans. Commun._ , vol. 69, no. 7, pp. 4732–4750, 2021.
  * [8] H. Sarieddeen, M.-S. Alouini, and T. Y. Al-Naffouri, “Terahertz-band ultra-massive spatial modulation MIMO,” _IEEE J. Sel. Areas Commun._ , vol. 37, no. 9, pp. 2040–2052, 2019.
  * [9] T. N. Do, G. Kaddoum, T. L. Nguyen, D. B. da Costa, and Z. J. Haas, “Multi-RIS-aided wireless systems: Statistical characterization and performance analysis,” _IEEE Trans. Commun._ , vol. 69, no. 12, pp. 8641–8658, Dec. 2021.
  * [10] C. Han and I. F. Akyildiz, “Distance-aware bandwidth-adaptive resource allocation for wireless systems in the terahertz band,” _IEEE Trans. Terahertz Sci. Technol._ , vol. 6, no. 4, pp. 541–553, 2016.
  * [11] Q. Wu, S. Zhang, B. Zheng, C. You, and R. Zhang, “Intelligent reflecting surface aided wireless communications: A tutorial,” _IEEE Trans. Commun._ , 2021.
  * [12] M. Di Renzo, K. Ntontin, J. Song, F. H. Danufane, X. Qian, F. Lazarakis, J. De Rosny, D.-T. Phan-Huy, O. Simeone, R. Zhang _et al._ , “Reconfigurable intelligent surfaces vs. relaying: Differences, similarities, and performance comparison,” _IEEE Open J. Commun. Soc._ , vol. 1, pp. 798–807, 2020.
  * [13] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 11, pp. 5394–5409, 2019.
  * [14] S. W. H. Shah, A. N. Mian, S. Mumtaz, A. Al-Dulaimi, I. Chih-Lin, and J. Crowcroft, “Statistical QoS analysis of reconfigurable intelligent surface-assisted D2D communication,” _IEEE Trans. Veh. Technol._ , vol. 71, no. 7, pp. 7343–7358, 2022.
  * [15] M. A. ElMossallamy, H. Zhang, L. Song, K. G. Seddik, Z. Han, and G. Y. Li, “Reconfigurable intelligent surfaces for wireless communications: Principles, challenges, and opportunities,” _IEEE Trans. Cognit. Commun. Networking_ , vol. 6, no. 3, pp. 990–1002, 2020.
  * [16] S. Dhok, P. Raut, P. K. Sharma, K. Singh, and C.-P. Li, “Non-linear energy harvesting in RIS-assisted URLLC networks for industry automation,” _IEEE Trans. Commun._ , vol. 69, no. 11, pp. 7761–7774, 2021.
  * [17] R. Hashemi, S. Ali, N. H. Mahmood, and M. Latva-aho, “Average rate and error probability analysis in short packet communications over RIS-aided URLLC systems,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 10, pp. 10 320–10 334, 2021.
  * [18] X. Li, Z. Xie, Z. Chu, V. G. Menon, S. Mumtaz, and J. Zhang, “Exploiting benefits of IRS in wireless powered NOMA networks,” _IEEE Trans. Green Commun. Networking_ , vol. 6, no. 1, pp. 175–186, 2022.
  * [19] J. Cheng, C. Shen, Z. Chen, and N. Pappas, “Robust beamforming design for IRS-aided URLLC in D2D networks,” _IEEE Trans. Commun._ , vol. 70, no. 9, pp. 6035–6049, 2022.
  * [20] F. Firyaguna, J. John, M. O. Khyam, D. Pesch, E. Armstrong, H. Claussen, H. V. Poor _et al._ , “Towards industry 5.0: Intelligent reflecting surface (IRS) in smart manufacturing,” _arXiv preprint arXiv:2201.02214_ , 2022\.
  * [21] M. Schellmann, “Capacity boosting by IRS deployment for industrial IoT communication in cm-and mm-wave bands,” in _2022 Joint European Conference on Networks and Communications & 6G Summit (EuCNC/6G Summit)_. IEEE, 2022, pp. 524–528.
  * [22] Y. Fang, S. Atapattu, H. Inaltekin, and J. Evans, “Optimum reconfigurable intelligent surface selection for wireless networks,” _IEEE Trans. Commun._ , 2022.
  * [23] C. Pan, H. Ren, K. Wang, J. F. Kolb, M. Elkashlan, M. Chen, M. Di Renzo, Y. Hao, J. Wang, A. L. Swindlehurst _et al._ , “Reconfigurable intelligent surfaces for 6G systems: Principles, applications, and research directions,” _IEEE Commun. Mag._ , vol. 59, no. 6, pp. 14–20, 2021.
  * [24] W. Tang, X. Chen, M. Z. Chen, J. Y. Dai, Y. Han, S. Jin, Q. Cheng, G. Y. Li, and T. J. Cui, “On channel reciprocity in reconfigurable intelligent surface assisted wireless networks,” _IEEE Wireless Commun._ , vol. 28, no. 6, pp. 94–101, 2021.
  * [25] E. Bjornson, O. Ozdogan, and E. G. Larsson, “Intelligent reflecting surface versus decode-and-forward: How large surfaces are needed to beat relaying?” _IEEE Wireless Commun. Lett._ , vol. 9, no. 2, pp. 244–248, Feb. 2020.
  * [26] Y. Liu, J. Zhao, M. Li, and Q. Wu, “Intelligent reflecting surface aided MISO uplink communication network: Feasibility and power minimization for perfect and imperfect CSI,” _IEEE Trans. Commun._ , vol. 69, no. 3, pp. 1975–1989, 2020.
  * [27] B. Al-Nahhas, A. Chaaban _et al._ , “Intelligent reflecting surface assisted MISO downlink: Channel estimation and asymptotic analysis,” in _GLOBECOM 2020_. IEEE, 2020, pp. 1–6.
  * [28] S. Hu, Z. Wei, Y. Cai, D. W. K. Ng, and J. Yuan, “Sum-rate maximization for multiuser MISO downlink systems with self-sustainable IRS,” in _GLOBECOM 2020_. IEEE, 2020, pp. 1–7.
  * [29] H. Shakhatreh, A. Sawalmeh, A. H. Alenezi, S. Abdel-Razeq, M. Almutiry, and A. Al-Fuqaha, “Mobile-IRS assisted next generation uav communication networks,” _arXiv preprint arXiv:2207.03622_ , 2022.
  * [30] L. Wang, Y. Chao, S. Cheng, and Z. Han, “An integrated affinity propagation and machine learning approach for interference management in drone base stations,” _IEEE Transactions on Cognitive Communications and Networking_ , vol. 6, no. 1, pp. 83–94, 2020.
  * [31] J. Cui, Y. Liu, Z. Ding, P. Fan, and A. Nallanathan, “Optimal user scheduling and power allocation for millimeter wave NOMA systems,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 3, pp. 1502–1517, 2017.
  * [32] M. N. Jung, C. Kirches, and S. Sager, “On perspective functions and vanishing constraints in mixed-integer nonlinear optimal control,” in _Facets of combinatorial optimization_. Springer, 2013, pp. 387–417.
  * [33] M. Schubert and H. Boche, “Iterative multiuser uplink and downlink beamforming under SINR constraints,” _IEEE Trans. Signal Process._ , vol. 53, no. 7, pp. 2324–2334, 2005.
  * [34] M. Lee and S. K. Oh, “A simplified iterative water-filling algorithm for per-user power allocation in multiuser MMSE-precoded MIMO systems,” in _VTC Spring 2008-IEEE Vehicular Technology Conference_. IEEE, 2008, pp. 744–748.
  * [35] Z. Wang, L. Liu, and S. Cui, “Channel estimation for intelligent reflecting surface assisted multiuser communications: Framework, algorithms, and analysis,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 10, pp. 6607–6620, 2020.
  * [36] S. Burer and A. N. Letchford, “Non-convex mixed-integer nonlinear programming: A survey,” _Surveys in Operations Research and Management Science_ , vol. 17, no. 2, pp. 97–106, 2012.
  * [37] D. Gale and L. S. Shapley, “College admissions and the stability of marriage,” _The American Mathematical Monthly_ , vol. 69, no. 1, pp. 9–15, 1962.
  * [38] Y. Xu and S. Mao, “User association in massive MIMO hetnets,” _IEEE Syst. J._ , vol. 11, no. 1, pp. 7–19, 2015.
  * [39] M. Najafi, V. Jamali, R. Schober, and H. V. Poor, “Physics-based modeling and scalable optimization of large intelligent reflecting surfaces,” _IEEE Trans. Commun._ , vol. 69, no. 4, pp. 2673–2691, 2020.
  * [40] Z. Mlika, E. Driouch, and W. Ajib, “User association under sinr constraints in HetNets: Upper bound and NP-hardness,” _IEEE Commun. Lett._ , vol. 22, no. 8, pp. 1672–1675, 2018.