A Pricing-Based Approach for Energy-Efficiency Maximization in RIS-Aided Multi-User MIMO SWIPT-Enabled Wireless Networks

In this work, we investigate the performance of a reconfigurable intelligent surface (RIS)-aided multi-user simultaneous wireless information and power transfer (SWIPT) network, where a multiple-input multiple-output (MIMO) base station (BS) serves multiple MIMO information receivers (IRs) while ensuring a minimum harvested power at multiple MIMO energy receivers (ERs). In order to improve the energy efficiency (EE) of the network, we consider a pricing-based performance metric called network utility. We then establish an optimization framework to jointly optimize the transmit precoding matrix (TPM) and phase shift matrix (PSM) to maximize the network utility function with constraints on the available transmit power at BS, minimum harvested power required at each ER, and unit modulus phase shift condition at RIS. Due to the non-convex nature of this problem, we divide it into two sub-problems where a sub-optimal solution of TPM and PSM are obtained separately using successive convex approximation and bisection search-based algorithms. Further, we propose an EE maximization (EEM) algorithm based on the block coordinate descent method to achieve the optimal solution of the master problem by iteratively obtaining the sub-optimal TPM, PSM, and network price using their respective algorithms. Moreover, we also prove that the solution obtained for each problem using their respective algorithm converges to the Karush-Kuhn-Tucker (KKT) optimum point of that problem. We also show the efficacy of the proposed algorithm using simulation results. In particular, we highlight the importance of using RIS in a multi-user MIMO SWIPT network and demonstrate the effect of various parameters on the network’s EE performance.


I. INTRODUCTION
The impending deployment of fifth-generation (5G) mobile communications around the world will be a major factor in driving productivity and will be the key enabler for long-envisaged verticals including personalised healthcare, manufacturing, smart energy grids, smart cities, finance, and transportation. However, realizing the ever-growing demand for a better communication network with improved quality of service (QoS) requirements such as lower power The associate editor coordinating the review of this manuscript and approving it for publication was Jie Gao . consumption, very high energy efficiency (EE), better spectral efficiency (SE), etc., researchers have already started to explore the evolution of 5G, commonly referred to as 5G and beyond (5GB) and sixth generation (6G) communications [1]. In this context, to enhance the experience of relaying in wireless communication networks, the idea of multiple passive reflecting surfaces/elements made of meta-materials has been floated to assist the communication between multiple devices such as base station (BS), users etc., [2]- [4]. This set of discrete reflecting elements is termed as reconfigurable intelligent surface (RIS) or intelligent reflecting surface (IRS), which is neither a part of the transmitter radio nor the receiver radio. It is a novel and low-cost/maintenance technique to control the wireless propagation medium, which until now had been deemed as uncontrollable. Accordingly, the fundamental role of a RIS is to affect the dispersion of wireless signals transmitted by other devices, without producing its signals [5]. The phase/reflection angle of each element in a RIS can be independently controlled or reconfigured through software. Due to their passive nature, RISs do not impose any thermal noise as they simply reflect the signals incident upon them. Therefore, RISs consume less power as compared to conventional decode-and-forward and amplify-and-forward relays, thus enhancing the EE of the communication network. Further, RIS uses smart reflection characteristics to create a virtual line of sight (LoS) link that helps in eliminating the fading caused by obstacles between the transmitter and receiver. In addition to this, RISs also provide multi-path propagation that leads to improvement in the rank of the channel along with an increase in the achievable diversity. Hence, a wireless network's performance can be significantly improved by mere adjustments of the angles of the reflective elements in a RIS [6], [7], which make it more remunerative than conventional relays for network operators.
While several key performance indices such as throughput, coverage, EE, and SE play pivotal roles in the design and deployment of RISs in a wireless network [8]- [12], this paper particularly focuses on the EE of the network 1 [13]. To achieve this, we not only consider the problem of maximizing the EE of the network but also implement simultaneous information and power transfer (SWIPT) [14], [15]. Accordingly, certain devices in the network have abilities to extract either only-power or only-symbol/information or some amount of power along with the desired symbol from the received signal using the energy harvesting (EH) technique. Now, both EE and EH are strictly related to the total power consumed by a network in the form of static and dynamic power to attain the desired level of QoS [16]. Static power consumption is the constant/fixed power required for non-communication-related tasks such as network signal processing, hardware maintenance, cooling, etc. On the other hand, dynamic power is the total transmit power used by the network to complete an end to end communication among several devices wirelessly. Due to continuous variations in wireless channel conditions, the total transmit power required by the network changes dynamically. Therefore, judicious selection of transmit power is required to significantly improve the network's performance.

A. RELATED LITERATURE AND MOTIVATION
A RIS-aided multiple-input single-output (MISO) network was considered in [10] and a two-timescale (TTS) transmission protocol was proposed to maximize the sum rate of the network. In particular, the transmit precoding vector was designed using instantaneous channel state information (CSI) and optimized phase shift matrix (PSM), which lead to a significant reduction in training overhead and design complexity. Similarly, in [11], joint optimization of transmitting precoder and PSM was performed to maximize the received power for a multiple RIS aided single user wireless network. Specifically, a semi-definite relaxation (SDR)-based algorithm was proposed to find a sub-optimal solution to the optimization problem. Next with regards to power efficiency, in [9], the authors minimized the total dynamic power by the joint design of active and passive beamformers for a RIS-aided MISO network. Similarly, joint symbol level precoding and PSM were optimized to minimize the total power usage in a RIS-aided multi-user network in [17]. In [18], joint optimization of the transmit power and PSM was performed to maximize the EE of a RIS-aided multi-user MISO network for Terahertz communication using covariance matrix adaptation evolution strategy and Dinkelbach's method.
While the works mentioned above are seminal for a RIS-aided network, the consideration of a MISO network or a network without EH capabilities or without an emphasis on EE, might not be ideal considering the requirements and specifications of future 6G networks. Accordingly, a multiuser MIMO network with multiple antennas at each node insinuates a more practical consideration. Similarly, the consideration of EH is also of paramount importance for future energy-efficient networks. Only a handful of works investigate networks with such considerations, where the communication between devices is assisted by RISs. For example, in [19] the trade-off between EE and SE was investigated in a RIS-aided multi-user MIMO uplink network and the resource efficiency was maximized by jointly optimizing the transmit precoding and RIS reflective beamforming. A similar RIS-aided multi-user MIMO uplink network was investigated in [20], where the passive beamforming and on-off reflecting modulation based information transfer were optimized to maximize the sum-rate. Next with regards to EH, the authors in [21] investigated the performance of SWIPT in a RIS-aided wireless multi-user network. Similarly, in [22], the authors investigated a RIS-assisted SWIPT multi-user network and maximized its sum rate by jointly optimizing the active and passive beamforming, albeit considering a single antenna at each user. However, to the best of the authors' knowledge no work to date has investigated the EE of a RIS-aided SWIPT enabled Multiuser MIMO Networks. Accordingly, the goal of this paper is to maximize the EE of a RIS-aided multi-user MIMO SWIPT network by jointly designing the optimal PSM and optimal precoding matrix.

B. CONTRIBUTIONS
In particular, we consider a multi-user MIMO SWIPT network having a multi-antenna BS serving multiple multi-antenna information receivers (IRs) while ensuring a minimum harvested power at multiple multi-antenna energy receivers (ERs). Unlike previous works, which have only considered sum rate maximization as a performance metric, we use a pricing-based approach and adopt a performance VOLUME 10, 2022 metric called network utility, which provides a striking balance between the sum rate and the total dissipated power, thus indirectly controlling the achievable EE of the network. We then provide a framework to jointly optimize the transmit precoding matrix (TPM) and PSM to maximize the network utility function. The primary distinctions of this work are summarized below.
• We formulate a network utility function maximization problem and jointly optimize TPM and PSM considering three different constraints: 1) maximum available transmit power available at BS, 2) minimum harvested power required at each ER, and 3) unit modulus phase shift at RIS. This problem is non-convex and extremely difficult to solve due to the involvement of phase and power constraints.
• We establish an optimization framework and reformulate the master problem into a simpler and more tractable form using the mean squared error (MSE) minimization approach. Then we use the block coordinate descent method to solve this reformulated problem by dividing it into two sub-problems and solve them separately.
In particular, first, we optimize TPM for fixed PSM and then optimize PSM for fixed TPM.
• To solve the TPM optimization problem, we propose a successive convex approximation (SCA) based algorithm. It uses the bisection search method to provide a near-optimal solution that converges at the Karush-Kuhn-Tucker (KKT) optimum point of this sub-problem. Similarly, we simplify the PSM optimization problem using the majorization-minimization method and propose a SCA based algorithm to find a near-optimal PSM that converges at the KKT optimum point.
• Further, we propose an EE maximization (EEM) algorithm based on the block coordinate descent method to achieve the optimal solution of the master problem by iteratively obtaining the sub-optimal TPM, PSM, and network price using their respective algorithms.
• Finally, the efficacy of the proposed algorithm is demonstrated using simulation results. In particular, we highlight the importance of using RIS in a multi-user MIMO SWIPT network and demonstrate the effect of various parameters on the network's EE performance.

1) ORGANIZATION
The remainder part of the paper is organized as follows: Section II discusses the considered RIS-aided multiuser wireless network in detail and formulated optimization problem. Sections III and IV present the problem formulation and its detailed solution, respectively. The simulation results are discussed in Section VI whereas Section VII concludes the work.

2) NOTATIONS
X and X * denote the converged solution and conjugate operator, respectively, of a given matrix X. tr(X), ||X|| F and |X| represent the trace, Frobenius norm and determinant, respectively. Re{x} denotes the real part of the complex value x. C a×b denote a complex vector of size a×b. The expectation operation is denoted by E{·}. (·) (·), (·) H and (·) T represent Hadamard product, Hermitian and transpose operations, respectively. arg{·} and diag(·) represents extraction of phase information and diagonalization operation, respectively. (·) † and (·) −1 denotes the matrix pseudo-inverse and inverse operations, respectively. I denotes the identity matrix. f x (x) is used for denoting the gradient of f w.r.t. the vector x and CN (0, σ 2 I) denote a random vector with zero mean and σ 2 variance.

II. SYSTEM MODEL
We consider a RIS-aided multi-user MIMO downlink wireless communication network which consists of a BS, N I information receivers (IRs) and N E energy receivers (ERs). We further assume that the BS is equipped with A B antennas, each IR and ER are fitted with A I and A E antennas, respectively. With the aid of N R reflective elements of RIS, BS transmits information to all IRs while providing sufficient energy to all ERs simultaneously using SWIPT protocol, as shown in the Fig. 1.

A. INFORMATION AND ENERGY TRANSFER
As mentioned earlier, BS is fitted with A B antennas to aid the signal transmission and it need to have channel state information (CSI) for applying transmit beamforming. Thus, for simplicity, we assume BS has perfect CSI using some standard estimation method 2 Therefore, using this CSI, BS designs an appropriate TPM and combines the symbols of all user to form a superimposed symbol given by where a n ∈ C c×1 represents the unit energy data symbol vector of the n th IR with c ∈ min(A B , A I ), A I represents the number of antennas at each IR. T n ∈ C A B ×c denotes the linear TPM corresponding to n th IR with ||T n || 2 F = P n and N I n=1 P n ≤ P max . P max is the maximum transmission power available at BS. Next, BS transmits this superimposed symbol to all users simultaneously. Thus, the signal received at the n th IR can be expressed as where = diag{φ 1 , φ 2 , . . . , φ N R } is the phase shift matrix (PSM) with φ i = e jθ i and θ i ∈ [0, 2π] denoting the phase shift of the i th reflective element of the RIS, L n,B ∈ C A I ×A B and L n,R ∈ C A I ×N R represent the channel gain from 2 CSI estimation in RIS-aided networks can be performed using some standard algorithms such as parallel factor decomposition (PARAFAC) [23]- [25]. However, similar to [26], [27], we assume perfect CSI in this paper for simplicity and analytical tractability. Moreover, on the reviewers suggestion and for better insights, we highlight the impact of CSI estimation error on the performance of the considered network in Fig. 11. BS and IRs, respectively, to the n th IR, P R,B ∈ C N R ×A B represents the channel gain for the link between BS and RIS, and d I ,n ∼ CN (0, σ 2 I I A I ) is the additive white Gaussian noise (AWGN). Substituting s from (1) into (2), b I ,n can be expressed as: whereL n = L n,B + L n,R P R,B . The corresponding achievable data rate of the n th IR is given by where K n is interference plus-noise covariance matrix given [28]. Further, the sum rate (SR) can be expressed as As mentioned earlier, this network also consists of low power ERs/sensor nodes which require certain amount of power (very low) to complete their designated task without any interruption. Therefore, considering this requirement, BS designs the beamforming/TPM so that alongwith satisfying the desired performance at each IRs, these ERs can also harvest sufficient power from received signal using SWIPT based EH [29]. Thus, using (1), the received signal at the j th ER can be expressed as where M j,B ∈ C A E ×A B and M j,R ∈ C A E ×N R are the channel gain from BS and RIS to the j th ER, respectively, and d E,j ∼ CN (0, σ 2 E I A E ) is the AWGN. Each ER applies EH protocol to harvest power from the received signal [14]. The total harvested power by the j th ER can be expressed as , N E } denotes the EH efficiency. 3 The weighted sum of the power harvested by all the ERs is given by H jM j , δ j is the energy weighting factor for the j th ER.

III. NETWORK UTILITY FUNCTION AND PROBLEM FORMULATION
Let us consider that BS consumes a fixed static power of P S B for performing non-communication task such as network signal processing, hardware maintenance, cooling, etc. Similarly, the static power consumed at the RIS is assumed to be P S I . Thus, total static power consumed by the considered network is given by Further, the total dynamic power consumed by the network is equal to the total transmit power used by BS. Therefore, from (1), the dynamic power can be evaluated as From (9) and (10), the total power consumed by the network is obtained as

A. NETWORK UTILITY FUNCTION
In this paper, we study the network's performance by investigating a pricing based EE maximization problem. Using (5) and (11), we define a new performance metric termed as network utility function (NUF) and expressed as where q ≥ 0 denotes the network price. Note that the power allocation problem will be equivalent to sum rate maximization problem when q → 0 and the power resource utilization cost becomes negligible. However, with increase in q, the optimal design of TPM and PSM becomes very important.

B. PROBLEM FORMULATION
As evident from (12), TPM T and PSM play a vital role in obtaining a desired performance. Therefore, it is extremely important to jointly optimize and design these parameters and maximize NUF. Thus, using (10), (8) and (12), we formulate a joint TPM and PSM optimization based NUF maximization problem considering the constraints of maximum available transmit at BS alongwith minimum power requirement of ERs and unit modulus of phase shift of RIS elements given by max T, whereH is the minimum power required by each ER. Clearly, due to coupling of variables T and , and involvement of EH constraint, this problem becomes non-convex in nature. Thus, it is extremely difficult to find its solution using standard methods.

IV. BLOCK COORDINATE DESCENT BASED ALGORITHM
As discussed above, problem (13) is extremely difficult to solve. Therefore, using minimum mean square error (MMSE) [33], problem (13) is simplified to a much simpler and tractable form. Thus, using (3) and approach similar to [33], the estimated signal vector at n th IR can be evaluated asâ where J n ∈ C A I ×c represents the decoding matrix. The mean square error (MSE) corresponding to (14) can be expressed as Using (4), (12) and (15) where V is the set of auxiliary matrices that denotes V = {V n ≥ 0, ∀n ∈ N I } and J is the set of J n for all IRs, and To find optimal J, we take first order partial derivative of l n (V, J, T, ) w.r.t. J n and equate it to zero. Thus, the optimal value of J n can be expressed as Similarly, the optimal value of V n can be expressed as After obtaining optimal decoding matrix J n and auxiliary matrix V n , we use these values to find sub-optimal TPM and PSM in next subsections. In particular, we use block coordinate descent (BCD) method to solve problem (16) and divide it into two sub problems where one problem will find sub-optimal TPM where as the other problem provides sub-optimal PSM.

A. SUB-OPTIMAL TPM
This section presents the optimization of TPM for fixed value of other variables (J, V and ). Using (15) and neglecting the constant terms, the maximization problem (16) can be reformulated in a minimization problem w.r.t T give as where Note that similar to previous problems, problem (20) is also non-convex. However, unlike previous problem, it is equivalent to a difference of convex (d.c.) program. Therefore, we can use successive convex approximation (SCA) method to solve problem (20) [34]. In particular, using first-order Taylor series expansion and Jensen' inequality, similar to [35], we obtain: where T (i) n is the previous value. After some algebra, we obtain 2Re tr Using (22), we reformulate the problem (20) as . Clearly, the OF in (23) is convex and it can be solved using standard optimization tool such as CVX [36]. However, the computational complexity of the CVX tool is very high. So, we further simply this problem and propose an algorithm using Lagrangian dual decomposition method that provides a near-optimal solution with much lower complexity [37]. Owing to the fact that it satisfies Slater's condition, we use the dual problem approach to find the solution as its dual gap is zero. The partial Lagrangian function of problem (23) can be expressed as where λ denotes the Lagrange multiplier. Next, using (24), we solve the following problem The Lagrangian function of problem (25) can be expressed as where ρ ≥ 0 denotes the dual variable. Differentiating L(T, ρ) w.r.t. T n and solving it for zero, we obtain Further, considering the complementary slackness condition of constraint (22), we obtain 2Re tr where λ i is the previous value of λ. Note that if inequality (28) hold true, then T n (λ i , 0) is the solution to the problem (25).
On the other hand, if inequality (28) does not satisfy, then the solution is given by T n (λ i , ρ o ) with Now, the dual problem is corresponding to g (λ) can be evaluated as If the inequality (31) holds true, then T n (0, ρ o ) is the solution to the problem (25). However, if inequality (31) does not satisfy, then the solution is given by T n (λ, ρ o ) [38], where λ is obtained by solving Due to involvement of ρ o in (29), obtaining a closed form solution of (32) is extremely difficult. So, in order to solve (32), we first obtain the behavior of P(λ) w.r.t. λ in the following lemma. Lemma 1: P(λ) decreases monotonically w.r.t. λ. Proof: Refer to Appendix A. Further considering the nature of P(λ), we propose Algorithm 1 that uses bisection search method to solve (23). Next using Algorithm 1 and SCA method, we propose Algorithm 2 to solve problem (20).
Theorem 1: Algorithm 2 provides a solution that converges at the KKT optimum point of (20).
Proof: Proof is similar to [39], so for paucity of space we have omitted the proof.  Evaluate P(λ) using (32) 12: if P(λ) ≥ P max then 13: λ in = λ. 14: else 15: λ fi = λ. Evaluate z(T (i) ) using (20) 5: Obtain {T n (i+1) , ∀n ∈ {1, 2 . . . , N I }} using Algorithm 1; 7: Evaluate z(T (i+1) ) using (23) 8: until i ≥ i max or |z(T (i+1) ) − z(T (i) )|/|z(T (i+1) )| ≤ ε where Y n = L H n,R J n V n J H n L n,R , U = P R,BT P H R,B , and W n = P R,BT H L H n,B J n V n J H n L n,R . C 1 is summation of entities that are independent of . Similarly, we obtain where B n = P R,B T n V n J H n L n,R and C 2 constitutes to the entities independent of . Further, using (8) with Y = N I n=1 Y n and X = N I n=1 W n − N I n=1 B n . Note that C 1 and C 2 are constant and have no impact on analysis, so we have neglected these terms. Using the identity from [40], we can write where H =H − Tr(M BT ) and ϒ = M R U T is a semidefinite matrix, as M R and U T are non-negative semidefinite matrices [40]. From (36) and (42), (37) is reformulated as where = Y U T . Similar to ϒ, is also a non-negative semidefinite matrix. Due to (42), problem (43) is a nonconvex problem. Therefore, similar to (20), we use SCA method to solve problem (43). Considering the fact that φϒφ is a convex function of φ, we have where φ (n) is the previous value. From (42) and (44), we obtain Thus using (45), the problem (43) can be transformed as whereĤ H +φ (n)H ϒφ (n) . We use the majorizationminimization algorithm to solve the problem (46) by subdividing it into simple and tractable sub-problems [41]. First we obtain the upper bound (ξ (φ|φ (n) )) of OF (χ (φ)) in (46) that satisfies the following conditions: We formulate a new sub-problem considering ξ (φ|φ (n) ) as the new OF which can be obtained using the following inequality [42] φ where Z = λ max I N R and λ max denotes the maximum eigen value of with . After some simplification, we obtain where Since, φ H φ = N R , thus φ H Zφ = N R λ max . Therefore, removing constant terms from (52), it is converted to a maximization problem given by where w (n) = (λ max I N R − )φ (n) − x * . Due to involvement of constraint (13d), problem (53) becomes a non-convex optimization problem. Also, unlike problem (23), its dual gap is not zero. Therefore, we use pricing based approach and reformulate the problem (53) as where p ≥ 0 denotes the price factor. The solution of this problem can be obtained as Considering the complementary slackness condition p K (p) −Ĥ = 0 of (45) with If (45) is satisfied, then p = 0 provide the optimal solution as φ(0). If (45) is not satisfied, then we need to find p for which the (55) gives the optimal solution. Therefore, using approach similar to (32), we propose Algorithm 3 that provides the optimal solution to problem (53).

Theorem 2: Algorithm 3 finds the optimal solution of problem (53) and problem (52).
Proof: For the proof refer to Appendix B. Using Algorithm 3, we propose Algorithm 4 that uses SCA method to solve the Problem (46). Theorem 3: Algorithm 4 provides a solution that converges at the KKT optimum point of (33).

C. OPTIMAL NETWORK PRICE (q ) FOR OPTIMAL EE
Using (5) and (11) The optimal price q that gives the maximum EE is obtained using following theorem. Theorem 4: The price q is the optimal price, if and only if the optimal precoder (T ) and optimal phase ( ) in problem (13) with respect to q satisfies the balance equation given by: Proof: Refer to Appendix D. Similar to [13], we use local maximizer of the problem (13) to obtain the local optimum value of q in the (i + 1) th iteration as where i is the previous iteration. Using (59), we find the optimal price q that can be obtained using iterative method. The proposed EEM algorithm is summarized in Algorithm 5.
D. SOLUTION TO PROBLEM (13) Using above discussion, we propose an iterative Algorithm 5 to solve (13). It uses Algorithm 2 and Algorithm 4 to find suboptimal T and , respectively, iteratively until convergence is achieved.
Proof: Refer to Appendix E. (5) and step (8), respectively. Therefore, these steps consumes the major portion of the computations with complexity discussed earlier in their respective section. Thus, the computational complexity of Algorithm 5 is given by

Algorithm 5 uses Algorithm 2 and Algorithm 4 in step
V. FEASIBILITY CHECK FOR PROBLEM (13) We formulate the following optimization problem to check the feasibility constraints (13b) and (13d) as follows: Problem (13) is feasible if and only if maximum OF of the above problem value is larger thanH . Due to coupling of TPM and PSM, obtaining an exact solution to the above problem is extremely difficult. Therefore, we obtain sub-optimal TPM and PSM by optimizing them alternately. Thus, first we obtain a sub-optimal TPM by solving following optimization problem The optimal solution to the above problem is obtained as . . , N I }, we have maximum value of OF in above problem as χP max , which represents the optimal energy beamforming [43].

VI. RESULTS AND DISCUSSION
In this section, the performance of the proposed algorithms and optimization framework are verified using exhaustive numerical simulations. Since ERs are energy constrained low power sensor nodes, thus, they are assumed to be placed in LoS of the BS and RIS. So, small scale fading of the link between RIS and ERs follows Rician distribution which is a combination of LoS and non-LoS (NLoS) components given byM and µ N R ϑ AoD where ϑ AoD riser,j ∈ [0, 2π] denotes the angle of departure and ϑ AoA riser,j ∈ [0, 2π] denotes the angle of arrival. γ and d represents the wavelength and antenna separation, respectively. Similar to (63), the link between BS to ERs and BS to RIS also follows Rician fading and are represented by M j,B ∈ C A E ×A B , ∀j ∈ {1, 2, . . . , N E } and P R,B ∈ C N R ×A B , respectively. Further, the link between BS to IRs and RIS to IRs are denoted by L n,B ∈ C A I ×A B and L n,R ∈ C A I ×N R , respectively, for all n ∈ {1, 2, . . . , N I }, and follow Rayleigh fading. Further, we consider a common large scale-path loss model for each link given by where PL 0 denotes the path loss at a reference distance D 0 and D j , ∀j ∈ {BSRIS, RISER, RISIR, BSIR, BSER} denotes the length for the links BS-RIS, RIS-ER, RIS-IR, BS-IR and BS-ER. δ j denotes the path loss exponent of the j th link. Also, (κ j , ∀j ∈ {1, 2, · · · , N E } = κ). Unless stated otherwise, the considered parameter values used for simulations are given in Table 1. Note that we consider without-RIS scenario (i.e., direct link only) as the benchmark to compare the proposed EEM algorithm. Fig. 2 represents the convergence behavior of the Algorithm 5 obtained for different number of RIS elements N R with P max = 10dBm. It can be observed that EE increases with increase in N R , which is evident from the fact that sum rate is directly proportional to N R (see (5)). Also, the obtained solution converges after a fixed number of iterations, which verifies Theorem 5. Also, it can be seen that the change in performance is around 150 − 200% for initial iterations. However, this change is around 0 − 2%, which is negligible compared to earlier change. This verifies the convergence of the algorithm graphically. Fig. 3 depicts the plot of average EE w.r.t. number of RIS elements N R obtained for three different cases: 1) Optimal PSM using Algorithm 5, 2) Random PSM, and 3) without RIS (wo-RIS, N R = 0). It can be seen that the performance is very low and constant when N R = 0. The main reason for this behavior is that with N R = 0, only direct link is used to decode information. Further, one can observe that there is a significant improvement in performance when N R increases. However, the achievable EE obtained using EEM algorithm is higher as compared to random PSM case. This highlights the importance of the resource optimization in the considered network in order to achieve maximum performance.
In Fig. 4, we plot the average EE versus the minimum harvested powerH for the three different cases with P max = 10dBm. As expected, the achievable EE reduces with increase inH . The reason for this behavior is that major portion of P max is used by ERs to achieve higherH . Therefore, net power received at each IR and hence the performance reduces with increase inH . It can also be seen that among three cases, EE is maximum for the optimal TPM and PSM case. For better understanding, we highlight the impact of energy harvesting efficiency (κ) on the EE of the considered network.   We plot average EE versus minimum harvested power at each ER (H ) in Fig. 5 for different values of harvesting efficiency at each user (K E = 4). In particular, we compare different/random efficiency (κ 1 = 0.3, κ 2 = 0.4, κ 3 = 0.5, κ 4 = 0.6) with equal efficiency (κ 1 = κ 2 = κ 3 = κ 4 = κ) at each user. As expected, performance increases with increase in efficiency because of increase in net harvested power.
All the above figures are obtained for fixed pathloss exponent of each link. However, depending on the environmental disturbances/obstacles such as buildings, trees, etc., the value of these parameters may change drastically. Therefore, to understand the impact of these pathloss exponent, we plot  average EE versus δ RIS in Fig. 6 considering δ RIS δ BSRIS = δ RISER = δ RISIR . It can be observed that because of the reduction in the strength of the reflected signal from RIS, the average EE for optimal and random case decreases with increase with δ RIS . Additionally, the performance is worst for the wo-RIS case as compared to other cases. However, for δ RIS ≈ 3, the performance of all the case are approximately same. Thus, using RIS in such scenarios has no significant impact on the network performance. Fig. 7 compares the average EE obtained after considering only direct link and only RIS links. Here, δ BSIR = δ BSER = δ D . As expected, the EE increases for each link with decrease in pathloss exponent. In addition to this, it can be seen that,  depending on δ D and δ RIS , the performance of RIS dominates the performance of direct link after a certain N R . For example, EE for RIS link is greater than that of the direct link when δ D = 3.6, δ RIS = 2 with N R > 100, and δ RIS = 2.1 with N R > 140. Similarly, EE for RIS link is greater than that of the direct link when δ D = 3.7, δ RIS = 2 ∀ N R , δ RIS = 2.1 ∀ N R , and δ RIS = 2.2 with N R > 100. Thus, however, increase in number of RIS elements can over come the double fading effect and provide better performance. Fig. 8 presents the variation of obtained EE w.r.t. location of ER (x ER ) with P max = 10dBm and BS at the origin. Similar to pathloss exponent, the performance for each case decreases with increase in distance between BS and ER. The main reason for this behavior is the reduction in effective power/strength of the received signal due to increase in distance. Further, there is a significant improvement in EE obtained using EEM algorithm as compared to other two cases. It can also seen that for high value of x ER , the performance of optimal and random cases are similar but still better than wo-RIS case. Therefore, in such case RIS provide a significant improvement over wo-RIS case.  In similar context, we show the impact of location of the IRs (x IR ) on the EE performance of all cases and plot average EE versus x IR in Fig. 9. Similar to Fig. 8, the performance for each case decreases with increase in distance between BS and IR, and optimal case provides the best performance. However, one can observe that the difference between the performance of optimal and random case is very low.
In Fig. 10, we study the impact of maximum transmit power (P max ) available at BS and observe the average EE gain for each scheme. It is obtained for three different values of number of elements at RIS N R . As expected, the average EE gain increases with increase in maximum transmit power for each N R . Also, N R = 50 provides much better performance than N R = 20 and N R = 10. Moreover, it can also be seen that after certain P max , there is no significant improvement in the achievable EE. Thus, one can choose P max carefully to avoid the wastage of power. Moreover, for better insights, we highlight the impact of CSI estimation error on the performance of the considered network in Fig. 11. It can be seen that average EE decreases due to increase in CSI estimation error.
Next, in Fig. 12, we plot the average EE versus number of antennas at BS (A B ) with P max = 10dBm and compare   the performance obtained using three cases discussed earlier.
It can be observed that with increase in A B there a significant improvement in achievable EE. This improvement is obtained because of the increase in diversity of system which increases linearly with increase in A B . In addition to this, one can also  observe that the difference between RIS-aided and wo-RIS cases is very high. Thus, using RIS in a multi antenna networks significantly improves the network performance.
Finally, we demonstrate the impact of the number of ERs (N E ) and IRs (N I ) on the network performance. Thus, we plot EE versus N E and EE versus N I in Fig. 13 and Fig. 14, respectively, for each case. It can be observed that the achievable EE increases with increase in the number of users in both the figures. However, observing both the graphs carefully for N E = N I , one can see that the EE obtained in Fig. 13 is lower than that of Fig. 14. The reason behind this difference is the minimum power requirement (EH limit) of ERs considered in this work. As the number of ER increases, their EH limit also increases. On the other hand, IR has no minimum power requirement.

VII. CONCLUSION
We investigated the performance of a RIS-adied multi-user MIMO SWIPT downlink network where a BS served information to multiple information receivers while ensuring a minimum harvested power at multiple energy receivers. A joint optimization of transmit precoding matrices and phase shift matrices was formulated for maximizing the network utility function considering constraints related to available transmit power available at BS, minimum harvested power required at each ER, and unit modulus phase shift condition at RIS. Due to non-convex nature of this problem, we divided it into two sub-problems and solve them separately. In particular, we optimized TPM and PSM separately using successive convex approximation and bisection search based algorithms. Further, we proposed an EEM algorithm based on block coordinate descent method that provided the optimal solution of the main problem by iteratively obtaining suboptimal TPM, PSM, and network price using their respective algorithms. Furthermore, we demonstrated the convergence behavior of the proposed algorithm using simulation results. Additionally, we also highlighted the impact of RIS in multiuser MIMO SWIPT networks. The importance of use of multiple antennas was also highlighted. Moreover, in order to avoid the wastage of power, the importance of judicious choice of available transmit power was also discussed.

APPENDIX D PROOF OF THEOREM 4
Sufficiency Part: Let Since T is the optimal precoder with respect to the optimal price q , it implies that q R T (n) , (n) P tot T (n) VOLUME 10, 2022 ≥ R T (n) , (n) P tot T (n) , ∀ T (n) , (n) ∈ S ∩ B, (82) where S depicts the feasible set of the problem (13) and B represents the norm ball set centered at T (n) , (n) with the radius r > 0. Thus, T (n) , (n) is the local maximizer for the priceq in the problem (13) because R T (n) , (n) −qP tot T (n) ≤ 0 = R T (n) , (n) −qP tot T (n) . (83) From (82) and (83), it can be concluded thatq is the optimal price such that the obtained solution T (n) , (n) can achieve the local maximum of the EE, and the balance equation is given as R T (n) , (n) −qP tot T (n) = 0. Necessity Part: Using the balance equation, we have R T (n) , (n) − q P tot T (n) ≤ R T (n) , (n) −q P tot T (n) = 0. (84) In other words, we can write R T (n) , (n)

APPENDIX E PROOF OF THEOREM 5
Firstly, the monotonic property of Algorithm 5 can be proved similar to [28]. Next, using KKT conditions corresponding to problem (20), its Lagrange function can be expressed as Here (92) corresponds to the chain rule and (95) is obtained after simplifying (19) using Woodbury matrix identity. From (95) and (90), we obtain ∇ T * n z(T, ) T n =T n = ∇ T * n R n (T, ) T n =T n .