SARIS: Scattering Aware Reconfigurable Intelligent Surface Model and Optimization for Complex Propagation Channels

The reconfigurable intelligent surface (RIS) is an emerging technology that changes how wireless networks are perceived, therefore its potential benefits and applications are currently under intense research and investigation. In this letter, we focus on electromagnetically consistent models for RISs inheriting from a recently proposed model based on mutually coupled loaded wire dipoles. While existing related research focuses on free-space wireless channels thereby ignoring interactions between RIS and scattering objects present in the propagation environment, we introduce an RIS-aided channel model that is applicable to more realistic scenarios, where the scattering objects are modeled as loaded wire dipoles. By adjusting the parameters of the wire dipoles, the properties of general natural and engineered material objects can be modeled. Based on this model, we introduce a provably convergent and efficient iterative algorithm that jointly optimizes the RIS and transmitter configurations to maximize the system sum-rate. Extensive numerical results show the net performance improvement provided by the proposed method compared with existing optimization algorithms.


I. INTRODUCTION
The reconfigurable intelligent surface (RIS) technology is gaining growing attention in academia and industry owing to its ability to turn the stochastic nature of the wireless propagation channel-always conceived as a black-box-into an optimizable variable [1].RISs are dynamic engineered metasurfaces, which are made of inexpensive scattering elements (unit cells) and low-complex/limited-power electronic circuits, which are spaced at sub-wavelength distances: the essentials are their dynamic reconfigurability in terms of scattering parameters performed in a nearly passive manner.In this context, adequate communication and channel models for RISs are of paramount importance to obtain a deep understanding of their achievable performance and to fully exploit their potential in future wireless networks [2].
Related work.In sub-wavelength implementations, there is an impelling need of considering the mutual coupling interactions among the unit cells: the discrete dipole approximation (DDA)-which is an analytical method for modeling the electromagnetic scattering from arbitrarily shaped objects [3], [4]-comes to help.In particular, the authors of [5] have recently introduced a mutual coupling-aware and electromagnetically consistent model for RIS with arbitrarily spaced unit cells, which are modeled as loaded wire dipoles.This model is based on the concept of mutually coupled antennas and accounts for near-and far-field channel conditions, thus offering great analytical tractability while ensuring electromagnetic consistency.Indeed, by capitalizing on the DDA, the polarizabilities of general natural or engineered materials (as is the case for RISs) are modeled by adjusting the parameters of the wire dipoles and the loads [3], [4], [6].
In [7], the authors have introduced a closed-form analytical formulation of the model in [5], while in [8] the model is exploited to derive the optimal RIS configuration that maximizes the received power in a single-user and singleantenna setup.A closed-form expression is given in the absence of mutual coupling, whereas a convergent iterative algorithm is proposed in the presence of mutual coupling.Such optimization framework is generalized in [9] for the multi-RIS multi-user multiple-input multiple-output (MIMO) interference channel.
With the exception of [9], such works are applicable to freespace propagation in the presence of an RIS.On the other hand, the authors of [9] model the multipath from environmental scattering objects (ESOs) by relying on a statistical fading model, which is superimposed (additive component) to the RIS-induced scattering.Inherently, this approach does not account for the interactions (e.g., multiple reflections) between the RIS and the scattering objects.Nevertheless, recent results show that multiple reflections between closelylocated scatterers (engineered or not) cannot be ignored [10].
Contributions.In this letter, we further extend the model in [5] to account for the multipath propagation originated from the presence of ESOs, which are modeled as a collection of loaded wire dipoles according to the DDA.In particular, the parameters of the dipoles modeling the ESOs (e.g., the length and the loads) are chosen to match the material properties of the objects [4], while the RIS loads are tunable on-demand.The proposed approach inherently captures the interactions between the RIS and the ESOs unlike the additive model in [9], which we prove is obtained as a particular case.
Based on the proposed model, we formulate a provably convergent optimization algorithm for computing the optimal load impedances of the RIS and the transmit precoding that maximize the sum mean squared error (SMSE) in a multipleinput single-output multi-user network, dubbed as SARIS.Compared with the state-of-the-art optimization algorithm introduced in [9], SARIS is shown to provide a higher sumrate while requiring fewer iterations and less execution time.
Notation.Matrices and vectors are denoted in bold font.(•) T , (•) H , and tr(•) stand for transposition, Hermitian transposition, and trace of a square matrix, respectively.• denotes the spectral norm and • F the Frobenius norm of a matrix.I N denotes the identity matrix of size N and 0 N ×M the all-zero N × M matrix, Lastly, j = √ −1 is the imaginary number.

II. SYSTEM MODEL
We consider a network comprising a transmitter with M antennas, L single-antenna receivers, and an RIS equipped with N reconfigurable (through tunable impedances) unit cells.Also, N s ESOs are available in the considered propagation environment.According to the DDA for natural (ESOs) or engineered (RIS) [3], [4] scatterers, the N unit cells and N s ESOs are modeled as loaded wire dipoles.The loads of the ESOs are not tunable, while the loads of the RIS are tunable on-demand for enhancing the communication performance.This allows us to adopt a unified model for RIS and ESOs.
Based on the model in [5], the end-to-end channel that accounts for the RIS and ESOs is formulated as [9, Eq. ( 6)] where {T, R, E} denote the transmitter, the receivers, and the scattering environment that includes the RIS and ESOs.
Compared with [5] and [9], the contribution of the ESOs is explicitly taken into account, as detailed next.Specifically, Z G ∈ C M ×M and Z L ∈ C L×L are the diagonal matrices containing the internal impedances of the voltage generators at the transmitter and the load impedances at the receivers, respectively; Z TT ∈ C M ×M and Z RR ∈ C L×L are the matrices containing the self and mutual impedances at the transmitter and receivers, respectively; and Z RT ∈ C L×M is the channel matrix of the direct link between the transmitter and receivers.The term Z RT − Z RE Z EE + Z SC −1 Z ET accounts for the signal scattered by the RIS and ESOs.Specifically, Z ET ∈ C (N +Ns)×M and Z RE ∈ C L×(N +Ns) are the channel matrices containing the mutual impedance between the transmitter and the scattering environment, and between the scattering environment and the receivers, respectively.Also, Z EE ∈ C (N +Ns)×(N +Ns) is the matrix containing the self and mutual impedances between the RIS and the ESOs, including the mutual coupling among the unit cells of the RIS.For further analysis, it is convenient to separate the contribution from the RIS and the ESOs.By denoting with the subscripts S and O the contributions of the RIS and the ESOs, respectively, the matrix Z EE can be partitioned into four blocks, as follows: By using a similar notation, the matrices Z ET and Z RE can be partitioned as follows: Lastly, the diagonal matrix Z SC contains the tunable load impedances of the RIS and the load impedances of the ESOs that model their material properties.Similar to Z EE , Z ET , and Z RE , the matrix Z SC can be partitioned as follows: where Z US ∈ C Ns×Ns is the diagonal matrix containing the loads of the ESOs and Z RIS ∈ C N ×N is the diagonal matrix containing the tunable impedances of the unit cells of the RIS.Each entry of the matrices of self and mutual impedances in (1) can be computed via the analytical frameworks in [7], [11], while the matrix Z US is determined by the material properties, e.g., the permittivity, of the ESOs, and is thus assumed to be fixed and given.On the other hand, the matrix Z RIS can be optimized, and, assuming that the RIS is nearly passive, is modeled as follows: where R 0 ≥ 0 is the resistance of each load, which accounts for the internal losses of the tuning circuits.If R 0 = 0, the RIS is lossless.The constraint R 0 ≥ 0 ensures that the RIS does not amplify the incident signal.It is customary to assume that R 0 is fixed, since it depends on the technology being used and the best choice is R 0 = 0 to minimize the losses.On the other hand, the reactance x i ∈ Q is tunable, with Q denoting the set (continuous or discrete) of feasible values for it.
To simplify the writing, we use the notation Hence, the received signal at the -th user equipment (UE), with = 1, . . ., L, is given by where W ∈ C M ×L denotes the precoding matrix at the transmitter, with W 2 F ≤ P and P being the power budget; s = [s 1 , . . ., s L ] T ∈ C L×1 is the transmit data vector, with E[ s 2 ] = 1; and n ∼ CN (0, σ 2 n I L ) is the noise vector at the receivers.Also, the following equivalent channels for the -th UE are introduced: with the lower-case bold letters denoting the -th row of the corresponding matrices in (1).With this notation, the sum-rate is given as follows: To simplify the notation and for further analysis, we introduce the following impedance matrices: As remarked in [8] and [9], the main difficulty when optimizing the sum-rate in (8) is given by the computation of the inverse matrix Z EE + Z SC −1 , since it is a non-diagonal matrix in the presence of mutual coupling.To deal with this matrix inversion, we invoke the Schur complement applied to block matrices [12].Specifically, assuming that Z OO is invertible, the matrix Z EE +Z SC −1 can be rewritten as shown in (15) at the top of this page.Thus, ( 1) and ( 7) can be simplified, respectively, as follows: The analytical expression of the end-to-end channel H E2E in ( 14) is conveniently formulated, since the optimization variables, i.e., the diagonal matrix Z RIS , is decoupled from the non-diagonal matrices Z SS and Z SOS that account for the mutual coupling between the unit cells of the RIS and for the interactions between the unit cells and the ESOs (including the mutual coupling among the wire dipoles of the ESOs), respectively.This facilitates the solution of optimization problems and provides deeper insights into the end-to-end channel.

A. Model Novelty: Modeling the ESOs
The analytical expression of H E2E in (14) clearly unveils the difference between the proposed model for the ESOs (i.e., the multipath) and closely related papers.Specifically, the authors of [9] have considered an additive statistical model for the multipath based on conventional fading models.The model in [9] can be retrieved as a special case from our proposed H E2E in (14) by setting Z SO = Z T OS = 0 N ×Ns , which yields Z ROS = −Z RS , Z SOS = 0 N ×N , and Z SOT = −Z ST .As a result, H E2E in (14) simplifies to The resulting end-to-end channel in ( 16) is given by the summation of the free-space channel in [5] and an additive multipath component given by Z RL Z RO Z −1 OO Z OT Z TG , which accounts for the scattering from the ESOs and is independent of the RIS, similar to the multipath model considered in [9].However, the channel in (16) inherently ignores the interactions between the RIS and the ESOs, thus resulting in less accurate modeling of the physical propragation environment (see, e.g., [10]) and suboptimal RIS designs that may fail to fully exploit its properties.Moreover, since the load impedances Z US can be arbitrarily chosen to mimic any natural material, while all the other self and mutual impedances depend mainly on the geometry of the scenario (e.g., see [7], [11]), the proposed model subsumes deterministic and statistical multipath channel models, i.e., in the case where the locations of the ESOs are assumed to be unknown (e.g., random).

III. PROBLEM FORMULATION AND SOLUTION
To optimize Z RIS subject to the constraint R 0 ≥ 0 in (5) and W ∈ C M ×L subject to the constraint W 2 F ≤ P , we consider the SMSE as objective function.The motivation for using the SMSE and its relation with the optimization of the sum-rate in (8) is detailed in [13].The SMSE is defined as Hence, the optimization problem is formulated as follows: (18) Due to the presence of the ESOs, the optimization problem in (18) is more challenging to solve as compared with those in [8] and [9].In fact, the matrix Z SOS is, in general, not a diagonal matrix even if the mutual coupling between the ESOs is negligible, i.e., Z OO is diagonal.This is because the matrices Z SO and Z OS are not diagonal matrices, in general.To tackle the problem in (18), we decouple it into two convex sub-problems that are solved iteratively, as detailed next.

A. Precoding Optimization
First, we solve the problem in (18) with respect to the optimization variable W while keeping Z RIS fixed.The resulting problem is convex and its solution is found by evaluating the KKT conditions.In this regard, let the Lagrangian and its gradient be respectively.Hence, the optimal closed-form solution is given by letting the expression in (19) to zero as [13, Eqs.(28), (29)]

B. RIS Optimization
Then, we solve the problem in (18) with respect to the optimization variable Z RIS while keeping W fixed.To this end, we devise an iterative algorithm that, similar to [8] and [9], leverages the repeated application of the Neumann series for efficiently tackling the inversion of Z SS + Z SOS + Z RIS .With respect to [8] and [9], however, we introduce an improved algorithm that accounts for the necessary conditions to make the Neumann series accurate by design.This is detailed next.Specifically, at the (i + 1)th iteration, the optimization variable Z RIS is set to is a diagonal matrix containing small improvements to the RIS configuration.With this approximation, we obtain where Therefore, we have As a result, the SMSE can be approximated as follows: The approximation in (21) is accurate provided that the condition ∆G i 1 is fulfilled at each iteration.In [8] and [9], this condition is taken into account by choosing very small values for the elements of δ.This approach, however, ignores the actual values of G i .To ensure that the Neumann series approximation is accurate by design at each iteration of the algorithm, we include the condition ∆G i 1 as a constraint of the optimization problem to be solved.This leads to the following simplified and robust problem formulation: where the constraint on |δ n | ensures that the Neumann series approximation is sufficiently accurate at each iteration.As in (19)-(20), the convex optimization problem in (26) is solved by setting the gradient of the Lagrangian to zero as  which has the following closed-form solution [13, Eq. ( 22)]: The normalization in (30) aims to strike a trade-off between the accuracy of the Neumann series approximation and the speed of convergence of the algorithm, i.e., the optimum is obtained when the inequality constraint in ( 26) is an equality.Hence, the small improvements δ are chosen as large as possible depending on the matrix G i at each iteration.
The complete algorithm that iterates between the solutions in ( 20) and (30), dubbed as SARIS, is provided in Algorithm 1. Specifically, once (30) is computed, the value of Z RIS is updated as Z i+1 RIS = Z i RIS + jIm{∆}, with ∆ = diag(δ H ), in order to preserve the constraint on the real part of Z RIS .Also, the imaginary part of Z i+1 RIS is projected onto the feasible set Q N .Compared with existing solutions, SARIS benefits from the inherent low-complexity of Algorithm 1, since the latter iterates between two simple closed-form expressions.

C. Computational Complexity and Convergence Analysis
The computational complexity of (20) and (30) is O(M 2 (M +3L)) and O(N (3N 2 +N (L+1)+3LM +M +L)), respectively.The complexity of (30) is dominated by the inversion of an N ×N matrix and is greater than the complexity of (20), since usually N M .Therefore, the total complexity of Algorithm 1 is O(N (3N 2 + N (L + 1) + 3LM + M + L)) multiplied by the number of required iterations to converge.
The proposed alternating optimization algorithm decouples the original non-convex problem in (18) into a series of convex sub-problems in the two optimization variables taken separately.Hence, the global objective function, i.e., the SMSE, is a non-increasing function across two consecutive iterations.Moreover, since it is, by definition, lower-bounded by zero, Algorithm 1 converges to a critical point of the problem in (18), which is a sufficient but not necessary condition for optimality [9], [13].

IV. NUMERICAL RESULTS
To evaluate the performance of Algorithm 1, we consider the scenario sketched in Fig. 1, where L = 2 UEs are present in the area of interest, unless otherwise stated.The simulation parameters are listed in Table I, where we remark that Z US = 0 models lossless metallic scattering objects.The inter-distance between the unit cells of the RIS is equal to d along both the xand y-axis.The number of unit cells along the xand y-axis is N x = N y = √ N , respectively.The locations of the transmitter, UEs, and the midpoint of the RIS are denoted by p BS , p UE1 , p UE2 , and p RIS , respectively.The N s ESOs are distributed in N c randomly located scattering clusters within a semicircle of radius R from the RIS.Each cluster comprises N O loaded wire dipoles, which are uniformly distributed (at random) within a disk of radius r centered at the locations of the clusters.The length of all wire dipoles (RIS and ESOs) is λ/2, where λ is the wavelength.All wire dipoles (RIS and ESOs) are aligned along the z-axis, i.e., Fig. 1 is a top view of the considered scenario, and all wire dipoles are orthogonal to the xy-plane (z = 0).The obtained sum-rate is averaged over 10 3 independent realizations for the locations of the scattering clusters and the loaded wire dipoles therein.
We compare the performance of SARIS against two benchmark schemes, namely the method in [9], which is denoted by BCD-wMMSE in the figures, and a mismatched approach.The latter is obtained by assuming as objective function H E2E in (16), which ignores the interactions between the RIS and the ESOs, and subsequenyly feeding the optimized W and Z RIS into the objective function in (18).In Fig. 2, we show the average sum-rate versus the number of iterations for different values of N and d, thus demonstrating the superior performance of SARIS both in terms of system sum rate and execution time to reach convergence, which is reported in Table .II.In Fig. 3, we illustrate the average sum-rate versus the number of scattering clusters for different values of N and d, demonstrating that ignoring the interactions between the RIS and the ESOs results in a significant degradation of the average sum-rate, especially as the number of ESOs increases.In Fig. 4 we show the average sum-rate versus the number of UEs L for different values of N and d.SARIS achieves higher performance for a small to moderate number of UEs, while the system becomes interference-limited for large L.
Finally, Fig. 5 shows the average sum-rate for small values of R 0 , assuming N = 64 and d = λ/8.In this case, the algorithm in [9] exhibits an unstable behavior.The reason is that the condition ∆G i 1 for the Neumann series to be accurate is not included in the formulation of the optimization problem, whereas it is taken into account at each iteration of Algorithm 1.This highlights the importance of adding the inequality constraint in (26), which makes the proposed approach robust-by-design when using the Neumann series approximation for small values of R 0 .For large values of R 0 , the two algorithms provide stable results, with the proposed algorithm converging faster.

TABLE I :
Simulation parameters.