Localization Performance Analysis and Algorithm Design of Reconfigurable Intelligent Surface-Assisted D2D Systems

Abstract

The research on high-precision and all-scenario localization using the millimeter-wave (mmWave) band is of great urgency. Due to the characteristics of mmWave, blockages make the localization task more complex. This paper proposes a cooperative localization system among user equipment (UEs) assisted by reconfigurable intelligent surfaces (RISs), which considers device-to-device (D2D) communication. RISs are used as anchor points, and position estimation is achieved through signal exchanges between UEs. Firstly, we establish a localization model based on this system and derive the UEs’ positioning error bound (PEB) as a performance metric. Then, a UE-RIS joint beamforming design is proposed to optimize channel state information (CSI) with the objective of achieving the minimum PEB. Finally, simulation analysis demonstrates the advantages of the proposed scheme over RIS-assisted base station positioning, achieving centimeter-level accuracy with a 10 dBm lower transmission power.

1. Introduction

In next-generation wireless networks, achieving high-precision localization in various scenarios remains a focus of further research. Current studies primarily rely on the accurate estimation of channel parameters, such as time of arrival (TOA), direction of arrival (DoA), or received signal strength indicator (RSSI), to enable localization by establishing communication links between transmitters and receivers. Among these, DOA-based methods such as multiple signal classification (MUSIC) [1] and estimation of signal parameters using rotational invariance techniques (ESPRIT) [2] algorithms, which utilize uniform linear arrays (ULA) for parameter estimation, have many applications. The positioning model based on RSSI has also been applied in recent research [3,4]. The RSSI value of radio waves during propagation is a function of the propagation distance, and this relationship can be used to estimate the distance. Utilizing RSSI based on pilot signal power for distance estimation at close range is a good choice, especially when there are no obstacles blocking the signal.

However, practical scenarios often involve obstacles between base stations (BSs) and user equipment (UEs) or mobile stations (MSs), resulting in non-line-of-sight (NLOS) situations. This becomes particularly prominent at millimeter-wave (mmWave) or higher frequencies, where signal penetration is weaker. Fortunately, with the advancements in metamaterials and radio frequency technology, reconfigurable intelligent surfaces (RISs) offer a promising solution to address NLOS localization [5]. Utilizing RISs for localization in mmWave scenarios has become a pivotal area of recent research. This is supported by studies such as [6,7,8], which demonstrate the capability of achieving centimeter-level positioning accuracy, even under conditions of obstructed links.

RIS is a revolutionary technology that intelligently reconfigures the wireless propagation environment by using a large number of low-cost passive reflecting elements integrated on a plane. It significantly enhances the performance of wireless communication networks [9]. RIS with multiple reflective elements provides higher beamforming (BF) gain and angular resolution, which is unmatched by ordinary scattering points. Despite the tremendous potential of RISs, from a communication design perspective, RIS faces new challenges in effectively integrating into wireless networks, such as reflection optimization, channel estimation, and deployment [9]. In recent years, researchers have made lots of meaningful work on RIS-aided communication positioning. The authors in [10] analyzed the impact of RIS-reflected paths on the accuracy of position estimation, but it only considered 2D scenarios. Other authors in [11,12] investigated the Cramér–Rao lower bound for absolute position estimation in 3D scenarios assisted by RIS, and optimized RIS-reflected BF to enhance localization performance. There are also considerations of cooperative localization methods. In [13,14], adding device-to-device (D2D) links between MSs was proposed to assist localization and the authors demonstrated that cooperative communication between MSs can provide additional localization gains through numerical results. However, the studies above are all based on the localization method with the BS as the necessary anchor point. Since RISs’ positions are known, it should be feasible to use RISs as anchor points for localization. In fact, the idea has been explored in [15,16]. The authors in [15] involved a single UE signal which is reflected back to itself through adjusting the RIS elements. After eliminating multipath effects, UE achieved localization itself. Ref. [16] introduced a D2D link based on [15], but its core is the same as [15]. In addition, the reflection design of the RIS is also a crucial aspect of positioning. Both [17,18] have conducted reflection BF designs by solely considering the RIS and utilizing discrete phase shifts. These studies have successfully demonstrated the effectiveness of such discrete phase shift designs. Hence, adopting a discrete design approach for localization is also a viable option.

Motivated by the above, this paper newly proposes a cooperative localization system among UEs with RISs as anchor points and conducts a performance analysis on this system. We consider both positioning relying solely on RIS and cooperative positioning among UEs. The proposed scheme differs significantly from previous papers in terms of details and introduces novel aspects in problem optimization. The main contributions of this paper are summarized as follows:

• We present a cooperative localization system utilizing RISs for millimeter-wave (mmWave) multi-input multi-output (MIMO) systems. In this system, UEs engage in D2D communication via RISs to facilitate localization, with the RISs serving as anchors without direct communication with BS or any access points. Leveraging RSSI and DOA methods such as ESPRIT, we extract parameters like angle and distance to establish geometric equations. The positions of UEs are then estimated by solving these equations. This proposed localization scheme not only significantly reduces the positioning pressure of the BS but also offers greater localization flexibility for UEs. Furthermore, it can serve as a supplementary means to enhance accuracy in BS-based UE localization.
• To evaluate the proposed system, we utilize the positioning error bound (PEB) of UEs as a performance metric. To accomplish this, we first calculate the Fisher information matrix (FIM) for the channel parameters. We then derive the transformation matrix between the FIMs of the channel parameters and the UEs’ positions. This allows us to determine the PEB for the UEs. With the aim of reaching the optimal performance of proposed system, we set the PEB as our objective. Since the received signal is significantly influenced by the channel state information (CSI), we innovatively propose a joint BF design of the transmitter UE’s BF and reflecting BF at RIS based on an alternating optimization algorithm. We illustrate our objective function and detail the algorithm steps of the proposed BF design to solve the function.
• The simulation results illustrate that, with the proposed BF design, our localization scheme achieves markedly superior positioning performance compared to scenarios without BF design or with random BF design. Moreover, by opting for smaller RIS sizes and an optimal number of quantization resolutions, our proposed scheme achieves positioning accuracy of ${10}^{-2}$ meters under the condition that signal-to-noise ratio (SNR) is 30 dB. This level of accuracy typically necessitates higher power consumption or larger RIS sizes in BS-based user localization systems.

The remainder of this paper is organized as follows. Section 2 introduces the system model, including the channel model, receiver model, and geometry relationship. Preliminary localization is presented, and the PEB is derived in Section 3. Section 4 demonstrates the optimization objective and proposes our algorithm. Section 5 shows the simulation analysis, including the impact of the antenna array size, the number of carriers, different quantization resolutions, and physical distance on PEB. There, we also conduct a comparative experiment. Conclusions are presented in Section 6.

Notations: Vectors are shown by bold lower-case letters and matrices by bold upper-case ones. $j=\sqrt{-1}$. $\mathbb{E}[·]$ denotes the expectation operator. $\Re \{·\}$ represents the real part of a complex variable. All vectors are column vectors by default. Transpose and Hermitian operations are denoted by ${(·)}^T$ and ${(·)}^H$, respectively. The Hadamard product is indicated by ⊙.

2. System Model

We consider a 2D localization scenario between UEs with two RISs assistance and no BS or access point, as presented in Figure 1, which can be called cooperative localization among UEs with RISs as anchors. It consists of RISs with L reflecting antenna elements in each ULA and two ULA users UE1 and UE2 with $N_{\mathrm{u}}$ antennas each. More users are feasible, but we aim to simplify the model and focus on the analysis of positioning performance. Two or more RIS units are necessary because we need a sufficient number of anchor points to construct the spatial topology of the devices. This conclusion is based on geometric intuition, and no formal proof is provided. Accordingly, we build the system model. The positions of UE1 and UE2 are, respectively, given by ${[q_x,q_y]}^T$ and ${[p_x,p_y]}^T$. As the system is operating in the mmWave band, near-field effects are neglected. The position of RISs, assumed to be points, are represented by ${[r{1}_x,r{1}_y]}^T$,${[r{2}_x,r{2}_y]}^T$. The rotation angles of the antenna arrays at UE1 and UE2 are denoted by ${\gamma }_q,{\gamma }_p\in [0,\pi ]$, respectively. The objective is to locate the positions of UE1 and UE2, given that the positions of RISs are fixed and known.

Consider the system where the working frequency of the carrier is $f_{\mathrm{c}}$ and the total bandwidth is B. In this system, each UE can control RISs and act as the transmitter for localization tasks. Thus, we assume that UE2 remains silent while UE1 transmits an Orthogonal Frequency Division Multiplexing (OFDM) signal with N subcarriers to UE2. Concurrently, UE2 receives both the direct signal from UE1 and the signal reflected from the RISs. During UE1’s transmission, UE1 remains silent as well.

2.1. Channel Model

The CSI matrix for the n-th subcarrier of the UE1-RIS1-UE2 route is formulated as

$${\mathbf{H}}_{\mathrm{RIS},1}\left[n\right]={\mathbf{H}}_{\mathrm{RR},1}\left[n\right]{\mathbf{\Omega }}_1{\mathbf{H}}_{\mathrm{TR},1}\left[n\right],$$

(1)

where ${\mathbf{H}}_{\mathrm{TR}}\left[n\right]$ is the CSI matrix of the UE1-RIS1 route, ${\mathbf{H}}_{\mathrm{RR}}\left[n\right]$ is the CSI matrix of the RIS1-UE2 route and $n\in \{-N/2,\cdots ,N/2\}$. $\mathbf{\Omega }=diag(e^{j{\omega }_1},e^{j{\omega }_2},\cdots ,e^{j{\omega }_L})\in {\mathbb{C}}^{L\times L}$ denote the reflection–coefficient matrix at the RIS, where ${\omega }_i\in [0,2\pi ),i=\{1,\cdots ,L\}$ are the reflection phase shifts. ${\mathbf{\Omega }}_1$ and ${\mathbf{\Omega }}_2$ denote the reflection–coefficient matrix at RIS1 and RIS2, respectively. ${\mathbf{H}}_{\mathrm{TR},1}\left[n\right]$ and ${\mathbf{H}}_{\mathrm{RR},1}\left[n\right]$ can be written as

$${\mathbf{H}}_{\mathrm{TR},1}\left[n\right]={\rho }_{11}{\mathit{\alpha }}_{\mathrm{in}}\left({\varphi }_{\mathrm{in},1}\right){\mathit{\alpha }}_{\mathrm{T}}^H\left({\alpha }_{\mathrm{t}1}\right)e^{-j2\pi B\frac{n}{N}{\tau }_{q,r1}},$$

(2)

and

$${\mathbf{H}}_{\mathrm{RR},1}\left[n\right]={\rho }_{12}{\mathit{\alpha }}_{\mathrm{R}}\left({\alpha }_{\mathrm{r}1}\right){\mathit{\alpha }}_{\mathrm{out}}^H\left({\varphi }_{\mathrm{out},1}\right)e^{-j2\pi B\frac{n}{N}{\tau }_{\mathrm{r}1,\mathrm{p}}},$$

(3)

where ${\rho }_{11}$ and ${\rho }_{12}$ are the free-space path loss occurred in the UE1-RIS1 and RIS1-UE2 link. ${\tau }_{\mathrm{q},\mathrm{r}1}$ and ${\tau }_{\mathrm{r}1,\mathrm{p}}$ are ToAs of the two links. ${\mathit{\alpha }}_{\mathrm{T}}\left({\alpha }_{\mathrm{t}1}\right)\in {\mathbb{C}}^{N_{\mathrm{u}}\times 1}$ and ${\mathit{\alpha }}_{\mathrm{R}}\left({\alpha }_{\mathrm{r}1}\right)\in {\mathbb{C}}^{N_{\mathrm{u}}\times 1}$ are the antenna array steering vectors and response vectors at UE1 and UE2, respectively. The i-th entry of ${\mathit{\alpha }}_{\mathrm{T}}\left({\alpha }_{\mathrm{t}1}\right)$ and ${\mathit{\alpha }}_{\mathrm{R}}\left({\alpha }_{\mathrm{r}1}\right)$ are ${\left[{\mathit{\alpha }}_{\mathrm{T}}\left({\alpha }_{\mathrm{t}1}\right)\right]}_i=e^{j2\pi (i-1)\frac{d}{\lambda }sin\left({\alpha }_{\mathrm{t}1}\right)}$ and ${\left[{\mathit{\alpha }}_{\mathrm{R}}\left({\alpha }_{\mathrm{r}}\right)\right]}_i=e^{j2\pi (i-1)\frac{d}{\lambda }sin\left({\alpha }_{\mathrm{r}1}\right)}$, respectively. ${\alpha }_{\mathrm{t}}$ is the AOD of UE1, and ${\alpha }_{\mathrm{r}}$ is the AOA of the UE2. d is the antenna element spacing of UE, and RIS1 is the same. $\lambda =c/f_c$, with c denoting the speed of light. ${\mathit{\alpha }}_{\mathrm{in}}\left({\varphi }_{\mathrm{in}}\right)$ and ${\mathit{\alpha }}_{\mathrm{out}}\left({\varphi }_{\mathrm{out}}\right)$ are the array response vectors of the UE1-RIS1 and RIS1-UE2 links, where ${\varphi }_{\mathrm{in}}$ is the AOA at the UE1-RIS1 link, and ${\varphi }_{\mathrm{out}}$ is the AOD at the RIS1-UE2 link. ${\mathit{\alpha }}_{\mathrm{in},1}\left({\varphi }_{\mathrm{in}}\right)$ is represented as ${\left[{\mathit{\alpha }}_{\mathrm{in}}\left({\varphi }_{\mathrm{in},1}\right)\right]}_i=e^{j2\pi (i-1)\frac{d}{\lambda }sin\left({\varphi }_{\mathrm{in},1}\right)}$, where $i\in \{0,1,\cdots ,L-1\}$, and ${\mathit{\alpha }}_{\mathrm{out}}\left({\varphi }_{\mathrm{out},1}\right)$ is defined similarly. Similar to (1), (2), and (3), the CSI matrix for the n-th subcarrier of the UE1-RIS2-UE2 route is formulated as

$${\mathbf{H}}_{\mathrm{RIS},2}\left[n\right]={\mathbf{H}}_{\mathrm{RR},2}\left[n\right]{\mathbf{\Omega }}_2{\mathbf{H}}_{\mathrm{TR},2}\left[n\right],$$

(4)

where

$${\mathbf{H}}_{\mathrm{TR},2}\left[n\right]={\rho }_{21}{\mathit{\alpha }}_{\mathrm{in}}\left({\varphi }_{\mathrm{in},2}\right){\mathit{\alpha }}_{\mathrm{T}}^H\left({\alpha }_{\mathrm{t}2}\right)e^{-j2\pi B\frac{n}{N}{\tau }_{\mathrm{q},\mathrm{r}2}},$$

(5)

and

$${\mathbf{H}}_{\mathrm{RR},2}\left[n\right]={\rho }_{22}{\mathit{\alpha }}_{\mathrm{R}}\left({\alpha }_{\mathrm{r}2}\right){\mathit{\alpha }}_{\mathrm{out}}^H\left({\varphi }_{\mathrm{out},2}\right)e^{-j2\pi B\frac{n}{N}{\tau }_{\mathrm{r}2,\mathrm{p}}}.$$

(6)

Moreover, UE2 will also send a positioning reference signal(PRS) to UE1 after UE2 performs this. The channel state information matrix of the UE2-RIS1-UE1 route is formulated as

$${\mathbf{G}}_{\mathrm{RIS},1}\left[n\right]={\mathbf{G}}_{\mathrm{RR},1}\left[n\right]{\mathbf{\Omega }}_1{\mathbf{G}}_{\mathrm{TR},1}\left[n\right],$$

(7)

where

$${\mathbf{G}}_{\mathrm{TR},1}\left[n\right]={\rho }_{12}{\mathit{\alpha }}_{\mathrm{out}}\left({\varphi }_{\mathrm{out},1}\right){\mathit{\alpha }}_{\mathrm{T}}^H\left({\alpha }_{\mathrm{r}1}\right)e^{-j2\pi B\frac{n}{N}{\tau }_{\mathrm{p},\mathrm{r}1}},$$

(8)

$${\mathbf{G}}_{\mathrm{RR},1}\left[n\right]={\rho }_{11}{\mathit{\alpha }}_{\mathrm{R}}\left({\alpha }_{\mathrm{t}1}\right){\mathit{\alpha }}_{\mathrm{in}}^H\left({\varphi }_{\mathrm{in},1}\right)e^{-j2\pi B\frac{n}{N}{\tau }_{\mathrm{r}1,\mathrm{q}}}.$$

(9)

${\mathit{\alpha }}_{\mathrm{T}}^H\left({\alpha }_{\mathrm{r}1}\right)$ and ${\mathit{\alpha }}_{\mathrm{R}}\left({\alpha }_{\mathrm{t}1}\right)$ are the antenna array steering vectors and response vectors at UE2 and UE1, respectively. They are defined similarly as ${\mathit{\alpha }}_{\mathrm{T}}\left({\alpha }_{\mathrm{t}1}\right)$ and ${\mathit{\alpha }}_{\mathrm{R}}\left({\alpha }_{\mathrm{r}1}\right)$. Similarly, the CSI matrix ${\mathbf{G}}_{\mathrm{RIS}2}\left[n\right]$ of the UE2-RIS2-UE1 route is expressed as

$${\mathbf{G}}_{\mathrm{RIS},2}\left[n\right]={\mathbf{G}}_{\mathrm{RR},2}\left[n\right]{\mathbf{\Omega }}_2{\mathbf{G}}_{\mathrm{TR},2}\left[n\right],$$

(10)

where

$${\mathbf{G}}_{\mathrm{TR},2}\left[n\right]={\rho }_{22}{\mathit{\alpha }}_{\mathrm{out}}\left({\varphi }_{\mathrm{out},2}\right){\mathit{\alpha }}_{\mathrm{T}}^H\left({\alpha }_{\mathrm{r}1}\right)e^{-j2\pi B\frac{n}{N}{\tau }_{\mathrm{p},\mathrm{r}1}},$$

(11)

$${\mathbf{G}}_{\mathrm{RR},2}\left[n\right]={\rho }_{21}{\mathit{\alpha }}_{\mathrm{R}}\left({\alpha }_{\mathrm{t}2}\right){\mathit{\alpha }}_{\mathrm{in}}^H\left({\varphi }_{\mathrm{in},2}\right)e^{-j2\pi B\frac{n}{N}{\tau }_{\mathrm{r}2,\mathrm{q}}}.$$

(12)

The CSI matrix of the UE1-UE2 route is formulated as

$${\mathbf{H}}_{\mathrm{UE}2}\left[n\right]={\rho }_3{\mathit{\alpha }}_{\mathrm{R}}\left({\theta }_{\mathrm{r}}\right){\mathit{\alpha }}_{\mathrm{T}}^H\left({\theta }_{\mathrm{t}}\right)e^{-j2\pi B\frac{n}{N}{\tau }_{\mathrm{u}}}.$$

(13)

The CSI matrix of the UE2-UE1 route is formulated as

$${\mathbf{G}}_{\mathrm{UE}1}\left[n\right]={\rho }_3{\mathit{\alpha }}_{\mathrm{T}}\left({\theta }_{\mathrm{t}}\right){\mathit{\alpha }}_{\mathrm{T}}^H\left({\theta }_{\mathrm{t}}\right)e^{-j2\pi B\frac{n}{N}{\tau }_{\mathrm{u}}}.$$

(14)

To obtain the distance of UE1 and UE2, each ULA needs to act as a single antenna (A1 and A2) to transmit signals. The CSI matrix of the A1-UE2 route is formulated as

$${\mathit{h}}_{\mathrm{UE}2}\left[n\right]={\rho }_3{\mathit{\alpha }}_{\mathrm{R}}\left({\theta }_{\mathrm{r}}\right)e^{-j2\pi B\frac{n}{N}{\tau }_{\mathrm{u}}}.$$

(15)

The CSI matrix of the A2-UE1 route is written as

$${\mathit{g}}_{\mathrm{UE}1}\left[n\right]={\rho }_3{\mathit{\alpha }}_{\mathrm{T}}\left({\theta }_{\mathrm{t}}\right)e^{-j2\pi B\frac{n}{N}{\tau }_{\mathrm{u}}},$$

(16)

where ${\rho }_3$ is the path loss in the UE1-UE2 link, and ${\tau }_{\mathrm{u}}$ is the path delay. ${\mathit{\alpha }}_{\mathrm{T}}\left({\theta }_{\mathrm{t}}\right)$ and ${\mathit{\alpha }}_{\mathrm{T}}\left({\theta }_{\mathrm{t}}\right)$ are defined in the similar way as ${\mathit{\alpha }}_{\mathrm{T}}\left({\alpha }_{\mathrm{t}1}\right)$.

2.2. Receiver Model

Assume that the transmitted PRS is $\mathbf{Fx}\left[n\right]\in {\mathbb{C}}^{N_{\mathrm{u}}\times 1}$ where $\mathbf{F}\in {\mathbb{C}}^{N_{\mathrm{u}}\times N_{\mathrm{u}}}$, $\mathit{x}\left[n\right]\in {\mathbb{C}}^{N_{\mathrm{u}}\times 1}$, ${\left|\right|\mathbf{Fx}\left[n\right]\left|\right|}_2=1$. $\mathbf{F}$ is transmit BF matrix. We assume that $\mathbf{F}$ at UE1 is denoted as ${\mathbf{F}}_1$, and $\mathbf{F}$ at UE2 is denoted as ${\mathbf{F}}_2$. According to (1)–(16), the signal ${\mathbf{y}}_1\left[n\right]$ reflected by the RISs and the signal ${\mathbf{y}}_3\left[n\right]$ directly transmitted from UE1 via a single antenna, both of which are received at UE2, can be expressed as

$$\begin{array}{c}\hfill {\mathbf{y}}_1\left[n\right]=\sqrt{P}({\mathbf{H}}_{\mathrm{RIS},1}\left[n\right]+{\mathbf{H}}_{\mathrm{RIS},2}\left[n\right]+{\mathbf{H}}_{\mathrm{UE}2}\left[n\right]){\mathbf{F}}_1\mathit{x}\left[n\right]+\mathbf{n}\left[n\right],\end{array}$$

(17)

$$\begin{array}{c}\hfill {\mathbf{y}}_3\left[n\right]=\sqrt{P}{\mathit{h}}_{\mathrm{UE}2}\left[n\right]+\mathbf{n}\left[n\right].\end{array}$$

(18)

The signal ${\mathbf{y}}_2\left[n\right],{\mathbf{y}}_4\left[n\right]\in {\mathbb{C}}^{N_{\mathrm{u}}\times 1}$ received at UE1 can be similarly expressed as

$$\begin{array}{c}\hfill {\mathbf{y}}_2\left[n\right]=\sqrt{P}({\mathbf{G}}_{\mathrm{RIS},1}\left[n\right]+{\mathbf{G}}_{\mathrm{RIS},2}\left[n\right]+{\mathbf{H}}_{\mathrm{UE}1}\left[n\right]){\mathbf{F}}_2\mathit{x}\left[n\right]+\mathbf{n}\left[n\right],\end{array}$$

(19)

$$\begin{array}{c}\hfill {\mathbf{y}}_4\left[n\right]=\sqrt{P}{\mathit{g}}_{\mathrm{UE}1}\left[n\right]+\mathbf{n}\left[n\right],\end{array}$$

(20)

where $\mathbf{n}\left[n\right]\in {\mathbb{C}}^{N_{\mathrm{u}}\times 1}$ is the additive white Gaussian noise drawn from $\mathcal{CN}(\mathbf{0},2{\sigma }^2\mathbf{I})$, and $\mathbf{I}$ is the identity matrix. P is the transmit power of the PRS.

2.3. Geometry Relationship

According to the geometric description of the model, the relationship between channel parameters and geometric information can be expressed as

$$\begin{array}{cc}\hfill {\tau }_{\mathrm{ri}}& ={\tau }_{\mathrm{q},\mathrm{ri}}+{\tau }_{\mathrm{ri},\mathrm{p}}=\frac{d_{\mathrm{ri},\mathrm{q}}}{c}+\frac{d_{\mathrm{ri},\mathrm{p}}}{c},\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill {\tau }_{\mathrm{u}}& =\frac{d_{\mathrm{q},\mathrm{p}}}{c},\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\theta }_{\mathrm{t}}& =\pi /2+arctan((p_x-q_x)/(p_y-q_y))+{\gamma }_q,\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill {\theta }_{\mathrm{r}}& =-\pi /2+arctan((p_x-q_x)/(p_y-q_y))+{\gamma }_p,\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {\alpha }_{\mathrm{ti}}& =arctan(({r\mathrm{i}}_y-q_y)/(q_x-{r\mathrm{i}}_x))+{\gamma }_q,\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill {\alpha }_{\mathrm{ri}}& =-arctan(({r\mathrm{i}}_y-q_y)/(q_x-{r\mathrm{i}}_x))+{\gamma }_p,\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill {\varphi }_{\mathrm{in},\mathrm{i}}& =-arctan(({r\mathrm{i}}_y-q_y)/(q_x-r{\mathrm{i}}_x)),\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\varphi }_{\mathrm{out},\mathrm{i}}& =arctan((p_y-{r\mathrm{i}}_y)/(p_x-r{\mathrm{i}}_x)),\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{i}}& ={\rho }_{\mathrm{i}1}{\rho }_{\mathrm{i}2}={\left(d_{\mathrm{ri},\mathrm{q}}\right)}^{-\mu /2}{\left(d_{\mathrm{ri},\mathrm{p}}\right)}^{-\mu /2},\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill {\rho }_3& ={\left(d_{\mathrm{q},\mathrm{p}}\right)}^{-\mu /2},\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill d_{\mathrm{ri},\mathrm{q}}& =\sqrt{{(q_x-r{\mathrm{i}}_x)}^2+{(q_y-r{\mathrm{i}}_y)}^2},\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill d_{\mathrm{ri},\mathrm{p}}& =\sqrt{{(p_x-r{\mathrm{i}}_x)}^2+{(p_y-r{\mathrm{i}}_y)}^2},\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill d_{\mathrm{p},\mathrm{q}}& =\sqrt{{(q_x-p_x)}^2+{(q_y-p_y)}^2},\hfill \end{array}$$

(33)

where $\mathrm{i}\in 1,2$, $\mu$ is the path loss exponent, and $d_{\mathrm{r},\mathrm{q}}$, $d_{\mathrm{r},\mathrm{p}}$, and $d_{\mathrm{p},\mathrm{q}}$ are the distances between the RISs and UE1, RISs and UE2, and UE1 and UE2, respectively.

3. Localization and Its Performance Metric

Based on the system model established in Section 2, we can estimate the parameters required for preliminary localization using the methods mentioned in the introduction. Subsequently, precision localization is achieved by adjusting the BF direction based on the preliminary localization results. We employ PEB as the evaluation metric for localization accuracy in this scheme. In Section 4, we will utilize this metric to formulate an objective function and propose a novel joint BF algorithm to optimize this objective function.

3.1. Preliminary Localization

Before optimizing the BF, we first roughly estimate the required parameters for localization. We employ the RSSI method for distance estimation and the ESPRIT method for receiving angle estimation. ${\mathbf{y}}_3\left[n\right]$ can also be written as:

$$\begin{array}{c}\hfill {\mathbf{y}}_3[n,l]=\sqrt{P}{\rho }_3{e}^{j2\pi l\frac{d}{\lambda }sin\left({\theta }_r\right)}{e}^{-j2\pi B\frac{n}{N}{\tau }_{\mathrm{u}}}+\mathbf{n}[n,l].\end{array}$$

(34)

Therefore,

$$\begin{array}{c}\hfill \widehat{{\rho }_3}=\frac{\frac{1}{N+1}{\sum }_{n=-N/2}^{N/2}\frac{1}{N_u}{\sum }_{l=0}^{N_u-1}\left|{\mathbf{y}}_3[n,l]\right|}{\sqrt{P}}.\end{array}$$

(35)

We set ${\mathbf{Y}}_1=[{\mathbf{y}}_1[-N/2],\cdots ,{\mathbf{y}}_1[N/2]]$ and calculate its covariance $\mathbf{R}=\frac{1}{N+1}{\mathbf{Y}}_1{\mathbf{Y}}_1^H$. Then, $\widehat{{\theta }_{\mathrm{r}}}$, ${\widehat{\alpha }}_{\mathrm{r}1}$, and ${\widehat{\alpha }}_{\mathrm{r}2}$ can be calculated by $\mathbf{R}$ with the aid of ESPRIT algorithm [2]. ${\widehat{\theta }}_{\mathrm{t}}$, ${\widehat{\alpha }}_{\mathrm{t}1}$, and ${\widehat{\alpha }}_{\mathrm{t}2}$ are calculated similarly. Based on the angles and distance, along with the known coordinates of the RISs, we can solve the equations in geometry relationship in Section 2. Thus, we depict the topological spatial location information between devices like Figure 1 and obtain the coordinates of the UEs.

3.2. Performance Metric for Localization

PEB is used to measure the localization scheme in a good number of studies. In this section, we derive PEB as a performance metric. As in [10,14], we use a two-stage method to obtain the FIM of the positions of UE1 and UE2. Firstly, we obtain the FIM of the channel parameters $\mathit{\eta }$. Then, the FIM for the UEs’ positions can be derived by computing a transformation matrix $\mathbf{T}$ that links the two FIM matrices. Finally, we use the FIM of the UEs’ positions to derive the Cramér–Rao Lower bound (CRLB) and position error bound for estimating the UE positions. According to the system model described above, all channel parameters $\mathit{\eta }$ can be stacked as

$$\mathit{\eta }={[{\rho }_3,{\theta }_{\mathrm{r}},{\alpha }_{\mathrm{r}1},{\alpha }_{\mathrm{r}2},{\theta }_{\mathrm{t}},{\alpha }_{\mathrm{t}1},{\alpha }_{\mathrm{t}2}]}^T,$$

(36)

We denote the unbiased estimator of $\mathit{\eta }$ as $\widehat{\mathit{\eta }}$. And firstly, we combine all the received signals into a new vector $y\left[n\right]$ defined as

$$\mathbf{y}\left[n\right]={[{\mathbf{y}}_1^T\left[n\right],{\mathbf{y}}_2^T\left[n\right],{\mathbf{y}}_3^T\left[n\right],{\mathbf{y}}_4^T\left[n\right]]}^T,$$

(37)

where $\mathbf{y}\left[n\right]\sim \mathcal{CN}(\mathit{\mu }\left[n\right],2{\sigma }^2\mathbf{I})$.

$$\mathit{\mu }\left[n\right]=[{\mathit{\mu }}_1^T\left[n\right],{\mathit{\mu }}_2^T\left[n\right],{\mathit{\mu }}_3^T\left[n\right],{\mathit{\mu }}_4^T\left[n\right]]$$

(38)

where ${\mathit{\mu }}_1\left[n\right]=\sqrt{P}({\mathbf{H}}_{\mathrm{RIS}1}\left[n\right]+{\mathbf{H}}_{\mathrm{RIS}2}\left[n\right]){\mathbf{F}}_1\mathit{x}\left[n\right]$, ${\mathit{\mu }}_2\left[n\right]=\sqrt{P}({\mathbf{G}}_{\mathrm{RIS}1}\left[n\right]+{\mathbf{G}}_{\mathrm{RIS}2}\left[n\right]){\mathbf{F}}_2\mathit{x}\left[n\right]$, ${\mathit{\mu }}_3\left[n\right]=\sqrt{P}{\mathbf{H}}_{\mathrm{UE}2}\left[n\right]{\mathbf{F}}_1\mathit{x}\left[n\right]$, ${\mathit{\mu }}_4\left[n\right]=\sqrt{P}{\mathbf{G}}_{\mathrm{UE}1}\left[n\right]{\mathbf{F}}_2\mathit{x}\left[n\right].$

Based on the Cramér–Rao theorem [19,20], the mean squared error of $\widehat{\mathit{\eta }}$ is bounded as

$$\mathbb{E}\left[(\widehat{\mathit{\eta }}-\mathit{\eta }){(\widehat{\mathit{\eta }}-\mathit{\eta })}^H\right]\ge {\mathbf{J}}_{\mathit{\eta }}^{-1},$$

(39)

where ${\mathbf{J}}_{\mathit{\eta }}\in {\mathbb{C}}^{12\times 12}$ is the FIM of $\mathit{\eta }$, and ${\left[{\mathbf{J}}_{\mathit{\eta }}^{-1}\right]}_{m,m}$ is the CRLB for the m-th parameter estimate. The $(m,n)$-th entry of ${\mathbf{J}}_{\mathit{\eta }}$ is defined as

$${\left[{\mathbf{J}}_{\mathit{\eta }}\left[n\right]\right]}_{m,n}=\mathbb{E}\left[-\frac{{\partial }^{2}lnp(\mathbf{y}\left[n\right];\mathit{\eta })}{\partial {\mathit{\eta }}_m\partial {\mathit{\eta }}_n}\right],$$

(40)

where $p\left(\mathbf{y}\right[n];\mathit{\eta })$ is the likelihood function of the random vector $\mathbf{y}\left[n\right]$ conditioned on $\mathit{\eta }$ and ${\mathit{\eta }}_m$ is the m-th entry of $\mathit{\eta }$. Since $\mathbf{y}\left[n\right]$ follows a symmetric complex normal distribution, the $(m,n)$-th entry for the FIM can be defined as

$${\left[{\mathbf{J}}_{\mathit{\eta }}\left[n\right]\right]}_{m,n}=\frac{P}{{\sigma }^2}\Re \left\{\frac{\partial {\mathit{\mu }}^H\left[n\right]}{\partial {\mathit{\eta }}_m}\frac{\partial \mathit{\mu }\left[n\right]}{\partial {\mathit{\eta }}_n}\right\},$$

(41)

The detailed calculation of $\mathit{\mu }\left[n\right]$ can be found in Appendix A.

Since our aim is to recover the UEs’ positions $\mathit{P}={[q_x,q_y,p_x,p_y]}^T$, we can obtain the FIM ${\mathbf{J}}_{\mathit{P}}$ associated with $\mathit{P}$ by properly transforming ${\mathbf{J}}_{\mathit{\eta }}$, which is shown as follows:

$${\mathbf{J}}_{\mathit{P}}\left[n\right]=\mathbf{T}{\mathbf{J}}_{\mathit{\eta }}\left[n\right]{\mathbf{T}}^T,$$

(42)

where the transformation matrix $\mathbf{T}\in {\mathbb{C}}^{4\times 12}$ is defined as:

$$\mathbf{T}=\frac{\partial {\mathit{\eta }}^T}{\partial \mathit{P}},$$

(43)

and the expression of ${\left[\mathbf{T}\right]}_i$ can be computed according to (21)–(33) and found in Appendix B.

Finally, we can obtain the entries of ${\mathbf{J}}_{\mathit{P}}$ as follows,

$${\left[{\mathbf{J}}_{\mathit{P}}\right]}_{m,n}=\sum _{n=-(N-1)/2}^{(N-1)/2}{\left[{\mathbf{J}}_{\mathit{P}}\left[n\right]\right]}_{m,n}=\sum _{n=-(N-1)/2}^{(N-1)/2}{\left[\mathbf{T}{\mathbf{J}}_{\mathit{\eta }}\left[n\right]{\mathbf{T}}^T\right]}_{m,n},$$

(44)

where $m,n\in \{0,1,\cdots ,4\}$. Therefore, the CRLB for estimating the UEs’ positions is the trace of the inverse matrix of ${\mathbf{J}}_{\mathit{P}}$: $\mathrm{CRLB}=\mathrm{tr}\left({\mathbf{J}}_{\mathit{P}}^{-1}\right)$. The PEB can be defined as

$$\mathrm{PEB}=\sqrt{\mathrm{CRLB}}=\sqrt{\mathrm{tr}\left({\mathbf{J}}_{\mathit{P}}^{-1}\right)}.$$

(45)

4. Optimization Objective and Solutions

According to the previous section, we know that the smaller PEB represents the better performance of the positioning system. Thus, in this section, our aim is to minimize the PEB. From Equations (44) and (45), it can be seen that PEB is mainly related to ${\mathbf{J}}_{\mathit{\eta }}$ and $\mathbf{T}$. During localization, the channel model and signal model are fixed, and $\mathbf{T}$ is a constant matrix. Therefore, PEB depends on the BF matrix and the reflection phase shifts at the RIS. Hence, we can minimize PEB by altering $\mathbf{F}$ and $\mathbf{\Omega }$.

Considering practical constraints, assume that the phase shift of each reflection element can only take discrete values. Let $Q_{\mathrm{r}}$ represent the number of control bits for each element, and suppose that the values of the discrete phase shifts are uniformly distributed in the range $(0,2\pi ]$. In $\mathbf{\Omega }$, ${\omega }_i\in \mathcal{F}$, where $\mathcal{F}=\{0,\frac{2\pi }{Q_r},\cdots ,\frac{2\pi (Q_{\mathrm{r}}-1)}{Q_{\mathrm{r}}}\}$. In addition, assume that the transmit BF matrix $\mathbf{F}=[f_1,\cdots ,f_{N_{\mathrm{u}}}]$ is constrained to take values from a codebook, where $f_k={[e^{j{v}_1},\cdots ,e^{j{v}_{N_{\mathrm{u}}}}]}^T$. We have $\mathbf{Fx}\left[n\right]=1/N_{\mathrm{u}}{[e^{j{v}_1},\cdots ,e^{j{v}_{N_{\mathrm{u}}}}]}^T$, where $v_j\in \mathcal{G}$, $\mathcal{G}=\{0,\frac{2\pi }{Q_{\mathrm{f}}},\cdots ,\frac{2\pi (Q_{\mathrm{f}}-1)}{Q_{\mathrm{f}}}\}$. Therefore, we can define an optimization objective function as

$$\begin{array}{cc}\hfill & \underset{w_i\forall i\in \{1,\dots ,L\},v_j\forall j\in \{1,\dots ,N_{\mathrm{u}}\}}{min}\sqrt{\mathrm{tr}\left({\mathbf{J}}_P^{-1}\right)}\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& {\omega }_i\in \mathcal{F},v_j\in \mathcal{G},\forall i=1,\dots ,L,\forall j=1,\dots ,N_{\mathrm{u}}.\hfill \end{array}$$

(46)

Since this objective function is non-convex, the block coordinate descent (BCD) algorithm can help us, which can be found in [14,21]. Hence, we propose a BCD-based joint alternative optimization as illustrated in Algorithm 1 and we call this optimization the joint BF design. In each iteration, we optimize ${\omega }_i$ while keeping ${\omega }_k$ fixed, where $k\in \{1,\cdots i-1,i+1,\cdots ,L^2\}$, and optimize $v_j$ while keeping $v_l$ fixed, where $l\in \{1,\cdots ,j-1,j+1,\cdots N_{\mathrm{u}}\}$. As a result, in ${\mathbf{\Omega }}_1$ or ${\mathbf{\Omega }}_2$, we have

$$\widehat{{\omega }_i}=arg\underset{w_i\in \mathcal{F}}{min}\sqrt{\mathrm{tr}\left({\mathbf{J}}_{\mathit{P}}^{-1}\right)},$$

(47)

where $i\in \{1,\dots ,L\}$, and in ${\mathbf{F}}_1$ or ${\mathbf{F}}_2$, we have

$$\widehat{v_j}=arg\underset{v_j\in \mathcal{G}}{min}\sqrt{\mathrm{tr}\left({\mathbf{J}}_{\mathit{P}}^{-1}\right)},$$

(48)

where $j\in \{1,\dots ,N_{\mathrm{u}}\}$. In order to approach the optimal solution as closely as possible, we alternate between optimizing $\mathbf{\Omega }$ while fixing $\mathbf{F}$ and optimizing $\mathbf{F}$ while fixing $\mathbf{\Omega }$. We calculate minimum PEB after optimizing $\mathbf{\Omega }$ and $\mathbf{F}$.

Finally, we present a concise complexity analysis of the proposed algorithm, which is pertinent given that many BCD-based algorithms undergo similar analyses. As delineated in Algorithm 1, during the i-th cycle of the t-th iteration, step 6 incurs a complexity of approximately $\mathcal{O}\left(N_u^3{L}^{4}N\right)$ operations. This complexity primarily stems from the multiplication and inversion of matrices. Notably, step 9 mirrors the operations of step 6. Consequently, the total complexity of the proposed algorithm is expressed as $\mathcal{O}\left(T(L{Q}_r+N_u{Q}_f)\left(N_u^3{L}^{4}N\right)\right)$, with T, typically less than ${10}^1$, representing the number of iterations for Algorithm 1.

Algorithm 1 The Proposed Joint Beamforming Design

1: $\mathrm{Objective}\mathrm{Function}:\mathrm{PEB}=\sqrt{\mathrm{tr}\left({\mathbf{J}}_{\mathit{P}}^{-1}\right)}$ 2: $\mathrm{Initialize}\mathit{\eta },{\mathbf{\Omega }}_0,{\mathbf{F}}_0,ϵ,{\mathrm{PEB}}_0\mathrm{and}\mathrm{set}t,i,j:=0.$ 3: while  $\Delta \mathrm{PEB}=|{\mathrm{PEB}}_t-{\mathrm{PEB}}_{t-1}|>ϵ$ do 4:     $t:=t+1$ 5:     for $i=1\mathrm{to}L$ do 6:         $\mathrm{Update}\widehat{{\omega }_i}\in \mathcal{F}\mathrm{s}.\mathrm{t}.\widehat{{\omega }_i}=argmin\sqrt{\mathrm{tr}\left({\mathbf{J}}_{\mathit{P}}^{-1}\right)}.$ 7:     end for 8:     for $j=1\mathrm{to}{N}_{\mathrm{u}}$ do 9:         $\mathrm{Update}\widehat{v_j}\in \mathcal{G}\mathrm{s}.\mathrm{t}.\widehat{v_j}=argmin\sqrt{\mathrm{tr}\left({\mathbf{J}}_{\mathit{P}}^{-1}\right)}.$ 10:     end for 11:     $\mathrm{Calculate}{\mathrm{PEB}}_t\mathrm{and}\Delta \mathrm{PEB}$ 12: end while 13: $\mathrm{Calculate}\mathrm{PEB}\mathrm{with}{\mathbf{\Omega }}_1^{*},{\mathbf{\Omega }}_2^{*}\mathrm{and}{\mathbf{F}}_1^{*},{\mathbf{F}}_2^{*}.$

5. Simulation Analysis

In this section, we analyze PEB to demonstrate the performance of the proposed system. Firstly, considering far field constraints, the position of the RISs are set at coordinates $[0,0]$ m and $[0,20]$ m, UE1 is at $[20,-10]$ m, and UE2 is at $[20,30]$ m. The signal-to-noise ratio (SNR) is defined as $\frac{P}{{\sigma }^2}$. If not specified, the default parameter settings are listed in

Table 1. Environmental parameters.

Parameters	Values	Description
L	32	The number of reflecting antenna elements for each RIS
$Q_{\mathrm{f}}$	6	The number of control bits for each UE
$Q_{\mathrm{r}}$	6	The number of control bits for each RIS
$N_{\mathrm{u}}$	32	The number of ULA antennas
P	0 dBm	The transmit power of PRS
${\gamma }_q$	$-\pi /9$	The rotation angle of UE1
${\gamma }_p$	$5\pi /12$	The rotation angle of UE2
N	31	The number of subcarriers
$ϵ$	${10}^{-6}$	The accuracy rating of Algorithm 1
d	0.005 m	The spacing between antenna arrays
$f_{\mathrm{c}}$	30 GHz	The center frequency of PRS
B	100 MHz	The subcarriers bandwidth
c	$3\times {10}^8\mathrm{m}/\mathrm{s}$	The velocity of light
$\mu$	2.08 [10]	The path loss exponent
${\sigma }^2$	−80 dBm [7]	The free space noise power

. These parameter settings have taken into account the avoidance of near-field effects.

5.1. Impact of the Joint BF Design

Firstly, we investigate the impact of the BF design on the transmitter and RIS on PEB. As shown in Figure 2, we compare four different cases, including (1) random $\mathbf{F}$ and random $\mathbf{\Omega }$; (2) random $\mathbf{F}$ and optimized $\mathbf{\Omega }$; (3) optimized $\mathbf{F}$ and random $\mathbf{\Omega }$; and (4) the joint optimized BF design based on the iterative solution of (46) by Algorithm 1.

The PEB performance of these four designs is depicted across various SNR levels in Figure 2. The results indicate that optimized $\mathbf{\Omega }$ or $\mathbf{F}$ yields superior PEB performance compared to random $\mathbf{\Omega }$ or $\mathbf{F}$, with PEB decreasing as SNR increases. This phenomenon arises from the ability to align transmitting beams, thereby concentrating them to achieve higher transmission gains. Similarly, without optimizing the phase design of the RIS, random reflection fails to attain higher acceptance gain for the receiving antenna. Of significant note, the joint BF design at the transmitter–RIS, based on Equation (44), achieves optimal PEB for the localization scheme, enabling centimeter-level localization accuracy under the given conditions.

5.2. Impact of $Q_r$ and $Q_f$

Next, we consider the impact of different phase quantization resolutions on the PEB generated by the joint BF design. As shown in Figure 3, it is evident that PEB decreases as $Q_{\mathrm{r}}$ and $Q_{\mathrm{f}}$ increase. However, it can also be observed that PEB growth saturates when $Q_{\mathrm{r}}\ge 6$ or $Q_{\mathrm{f}}\ge 6$.

5.3. Impact of L, $N_u$, and N

Next, we consider the impact of subcarrier count design and the size design of RIS on the localization performance. As shown in Figure 4 and Figure 5, it is obvious that with the increase of L, $N_{\mathrm{u}}$, and N, a significant reduction in PEB can be observed.

This is intuitive, as a larger RIS represents more RIS reflecting elements, thereby providing more reflective links to deliver the required localization information. Having longer antennas on the UEs implies the capability to transmit and receive more link information, naturally aiding in improving localization accuracy. Similarly, a larger N implies more subcarriers, which also brings more links capable of providing measurement information. The results further indicate that by selecting an appropriately sized RIS, achieving centimeter-level localization accuracy can be easily realized with a large number of reflecting elements.

5.4. Impact of the Distance between UEs

Furthermore, the distance between UE1 and UE2 should also be considered. Constrained by the signal attenuation in free space, PEB gradually increases as the distance between devices grows larger in Figure 6. It is also noted that increasing power can partially compensate for this loss of accuracy.

5.5. Comparison

Finally, we compare the performance gap between the proposed scheme and the scheme (Comparison) that BS locates UEs with the assistance of RISs as in Figure 7. A similar model can be seen in [11] for the localization scheme used for comparison. It is essential to note that the comparative scheme under consideration lacks D2D links. Furthermore, in practical scenarios, the transmission power of the BS, at $P_{\mathrm{t}}=10\mathrm{dBm}$, is significantly higher than that of the UEs. Typically, the distance $d^{\prime }$ between the BS and the UEs should far exceed the distance between the UEs themselves. Accordingly, we have set $d^{\prime }=5\times d_{\mathrm{p},\mathrm{q}}$. The BS is positioned at coordinates [20, −170] meters and is equipped with a single antenna, hindered by obstacles. Localization can only be achieved through reflection with the assistance of RISs. Other parameters were set to be as consistent as possible.

As can be seen from Figure 8, the proposed scheme clearly outperforms the traditional RIS-assisted localization and lower power can yield improved performance. Additionally, replacing RISs with a larger size or adding more RISs can bridge the performance gap between the schemes. This demonstrates the feasibility and advantages of cooperative localization schemes. However, it should be noted that our proposed scheme requires collaborative transmission of signals among devices, whereas the comparative scheme only requires a single transmission. As a result, our scheme incurs twice the complexity. Moreover, this D2D localization scheme distributes the localization burden to UEs rather than BS, reducing the power pressure on BS. Even as a performance supplement for conventional localization schemes, it remains a viable option, as can be seen in

Table 2. Difference between our proposed and comparison.

Points	Our Proposed	Comparison
Distance between transmitter and receiver	40 m	200 m
Transmission power	0 dBm	10 dBm
Transmission path	not blocked and collaborative	blocked
Localization performance	better	good
Time complexity	double the comparison	normal

.

6. Conclusions

In this paper, we proposed a system for D2D cooperative localization in mmWave MIMO systems using RISs as anchor points, without relying on BS. UEs transmit PRS to each other through RISs to achieve localization. A system model is then established based on this scheme. The PEB, which evaluates the coordinate error of the two UEs’ positions, is derived, and a joint transmitter–RIS BF design algorithm based on alternating optimization is proposed to minimize the PEB. The simulation results demonstrate that the proposed algorithm significantly improves the localization performance and the proposed scheme can achieve centimeter-level localization accuracy.

In conclusion, our proposed cooperative localization scheme leveraging RISs presents a promising solution for enhancing localization accuracy and flexibility in mmWave MIMO systems. Future research directions may include further optimization algorithms and practical implementations to validate the effectiveness of our proposed scheme in real-world scenarios.