Gui Zhou
email: g.zhou@qmul.ac.uk

Cunhua Pan
email: c.pan@qmul.ac.uk

Hong Ren
email: hren@seu.edu.cn

Kezhi Wang
email: kezhi.wang@northumbria.ac.uk

Fellow, IEEE Marco Di Renzo
email: marco.di-renzo@universite-paris-saclay.fr

) G Zhou

Fairness-Oriented Multiple RIS-Aided mmWave Transmission: Stochastic Optimization Methods

Keywords: Reconfigurable intelligent surface (RIS), intelligent reflecting surface (IRS), millimeter wave communications, stochastic optimization, robust beamforming design

In millimeter wave (mmWave) systems, it is challenging to ensure the reliable connectivity of communication due to the high sensitivity to the presence of blockages. In order to improve the robustness of mmWave systems under the presence of random blockages, multiple reconfigurable intelligent surfaces (RISs) are deployed to enhance the spatial diversity gain, and robust beamforming is then designed based on a stochastic optimization for minimizing the maximum outage probability among multiple users to ensure fairness. Under the stochastic optimization framework, we adopt the stochastic majorization-minimization (SMM) method and the stochastic successive convex approximation (SSCA) method to construct deterministic surrogate problems at each iteration for new channel realizations, and to obtain closed-form solutions of the precoding matrix at the base station (BS) and the passive beamforming vectors at the RISs. Both stochastic optimization methods are proved to converge to the set of stationary points of the original stochastic problems. Finally, simulation results show that the proposed robust beamforming in RIS-aided systems can effectively compensate for the performance loss caused by the presence of random blockages, especially when the blockage probability is high, compared with benchmark solutions.

I. INTRODUCTION

Millimeter wave (mmWave) communication is expected to be a promising technology to meet the growing demand for data rate in current and future wireless networks due to the abundant spectrum available at high frequencies. High frequency communication inevitably has severe attenuation, but this can be compensated by an array of reasonable size containing a large number of antennas due to the small wavelength [START_REF] Rappaport | Millimeter wave mobile communications for 5G cellular: it will work![END_REF]. In addition, the high-directional beams of large arrays are capable of mitigating the inter-user interference. However, mmWave communication systems suffer from high penetration loss [START_REF] Raghavan | Spatio-temporal impact of hand and body blockage for millimeter-wave user equipment design at 28 GHz[END_REF]- [START_REF] Raghavan | Statistical blockage modeling and robustness of beamforming in millimeter-wave systems[END_REF]. Hence, mmWave systems are much more susceptible to the presence of spatial blockages than sub-6 GHz systems, and the reliability of the communication links for the entire network cannot always be guaranteed [START_REF] Raghavan | Spatio-temporal impact of hand and body blockage for millimeter-wave user equipment design at 28 GHz[END_REF]- [START_REF] Raghavan | Statistical blockage modeling and robustness of beamforming in millimeter-wave systems[END_REF].

In particular, spatial blockages can be divided into static blockages (e.g., buildings and other immobile fixtures), dynamic blockages (e.g., humans, vehicles, or moving obstructions) and self-blockages (e.g., hand blocking of the user itself and blockage from other body parts). Some statistical models are established to characterize the properties of the random dynamic blockages and self-blockages [START_REF] V. Raghavan | Antenna placement and performance tradeoffs with hand blockage in millimeter wave systems[END_REF]- [START_REF] Bai | Analysis of blockage effects on urban cellular networks[END_REF]. The authors of [START_REF] Akdeniz | Millimeter wave channel modeling and cellular capacity evaluation[END_REF] have developed a distance-dependent blockage probability model, in which the probability of link blocking increases exponentially with the link length. Furthermore, [START_REF] Nishio | Proactive received power prediction using machine learning and depth images for mmWave networks[END_REF] and [START_REF] Alrabeiah | Deep learning for mmWave beam and blockage prediction using sub-6 GHz channels[END_REF] proposed a mechanism to predict the blockage probability via machine learning. When the blockage probability is obtained, some robust beamforming design strategies have been proposed in the recent literature [START_REF] Kumar | Rate and reliability trade-Off for mmWave communication via multi-point connectivity[END_REF], [START_REF] Iimori | Stochastic learning robust beamforming for millimeter-wave systems with path blockage[END_REF] to address channel uncertainties caused by the presence of random blockages. Specifically, the authors of [START_REF] Kumar | Rate and reliability trade-Off for mmWave communication via multi-point connectivity[END_REF] proposed a worstcase robust beamforming design for application to coordinated multipoint (CoMP) systems in which all possible combinations of blockage patterns were considered. Due to the high computational complexity of the ergodic method in [START_REF] Kumar | Rate and reliability trade-Off for mmWave communication via multi-point connectivity[END_REF], an outage-minimum strategy based on a stochastic optimization method was proposed for CoMP systems in [START_REF] Iimori | Stochastic learning robust beamforming for millimeter-wave systems with path blockage[END_REF] to improve the robustness of mmWave systems. The idea of adopting multiple base stations (BSs) in CoMP systems is effective to compensate for the performance loss caused by the presence of random blockages by exploiting spatial diversity gains. However, this will incur excessive hardware cost and power consumption. Another promising scheme proposed in [START_REF] Zhou | Stochastic learning-based robust beamforming design for RIS-aided millimeter-wave systems in the presence of random blockages[END_REF] is to deploy cost-efficient reconfigurable intelligent surfaces (RISs) in mmWave systems to create an alternative communication link via the RISs.

Due to the promising advantages in terms of providing energy-and spectral-efficient communications, RISs have attracted extensive research attention from both academia and industry [START_REF] Renzo | Smart radio environments empowered by reconfigurable AI meta-surfaces: An idea whose time has come[END_REF]- [START_REF] Yu | IRS-assisted green communication systems: Provable convergence and robust optimization[END_REF]. An RIS is a thin surface consisting of nearly-passive and reconfigurable reflecting elements, which reflects the impinging radio waves without adopting radio frequency (RF) chains.

The passive elements on the RIS can be tuned to alter their electromagnetic response such that the signals reflected from an RIS can be constructively superimposed to enhance the signal power at the intended receiver or destructively combined to avoid the information leakage to undesired receivers. These characteristics make RISs appealing to be applied in various communication systems as shown in [START_REF] Pan | Reconfigurable intelligent surface for 6G and beyond: Motivations, principles, applications, and research directions[END_REF]. For instance, RISs can be applied in single-cell multiple-input and multiple-output (MIMO) systems [START_REF] Pan | Intelligent reflecting surface aided MIMO broadcasting for simultaneous wireless information and power transfer[END_REF]- [START_REF] Zhi | Power scaling law analysis and phase shift optimization of RIS-aided massive MIMO systems with statistical CSI[END_REF], multicell MIMO communications [START_REF] Pan | Multicell MIMO communications relying on intelligent reflecting surfaces[END_REF], simultaneous wireless information and power transfer (SWIPT) systems [START_REF] Pan | Intelligent reflecting surface aided MIMO broadcasting for simultaneous wireless information and power transfer[END_REF], [START_REF] Wu | Joint active and passive beamforming optimization for intelligent reflecting surface assisted SWIPT under QoS constraints[END_REF], secure communications [START_REF] Hong | Artificial-noise-aided secure mimo wireless communications via intelligent reflecting surface[END_REF], mmWave systems [START_REF] Mursia | RISMA: Reconfigurable intelligent surfaces enabling beamforming for IoT massive access[END_REF]- [START_REF] Perović | Channel capacity optimization using reconfigurable intelligent surfaces in indoor mmWave environments[END_REF] and THz systems [START_REF] Pan | Sum rate maximization for intelligent reflecting surface assisted terahertz communications[END_REF], [START_REF] Ning | Channel estimation and transmission for intelligent reflecting surface assisted THz communications[END_REF].

Although the performance advantages of deploying an RIS in mmWave systems have been demonstrated in recent contributions, there still exist major open problems to solve. The authors of [START_REF] Di | Practical hybrid beamforming with finite-resolution phase shifters for reconfigurable intelligent surface based multi-user communications[END_REF] only considered the BS-RIS-user channels and assumed that the direct BS-user communication links were completely blocked by obstacles. However, this assumption only applies to the case of static blockages with outage probability of 1, but not to the case of dynamic blockage with outage probability varying between 0 and 1 [START_REF] Raghavan | Spatio-temporal impact of hand and body blockage for millimeter-wave user equipment design at 28 GHz[END_REF]- [START_REF] Raghavan | Statistical blockage modeling and robustness of beamforming in millimeter-wave systems[END_REF]. The numerical results in [START_REF] Wang | Intelligent reflecting surface-assisted millimeter wave communications: Joint active and passive precoding design[END_REF] showed that the gain from additional reflection channels could compensate for the performance loss caused by the presence of random blockages, but the impact of blockages was not considered in the beamforming design. Most recently, we have considered the robust beamforming design for RIS-aided mmWave communication systems in [START_REF] Zhou | Stochastic learning-based robust beamforming design for RIS-aided millimeter-wave systems in the presence of random blockages[END_REF] by taking the random blockages into consideration. However, the objective function therein is to minimize the sum outage probability, which cannot ensure the fairness for all the users.

A. Novelty and contributions

Against the above background, this paper proposes a robust transmission strategy in an RISaided mmWave communication system to deal with the channel uncertainties caused by the random blockages while ensuring the fairness among the users. Some common approaches to handle the imperfect and partial CSI problem are the outage constrained robust optimization and the worst-case robust optimization techniques [START_REF] Zhou | A framework of robust transmission design for IRS-aided MISO communications with imperfect cascaded channels[END_REF]. However, both of them still need the estimation of the instantaneous CSI, and the worst-case robust method is conservative and hence suboptimal due to the low probability of occurrence for the worst case. An alternative approach is to design the robust beamforming by optimizing the statistical performance under a stochastic optimization framework in which only the long-term CSI is required. In what follows, we propose a maximum outage probability minimization problem to be solved by using the stochastic optimization framework. Specifically, the main contributions of this work are summarized as follows:

• To the best of our knowledge, this is the first work exploring the robust beamforming design for an RIS-aided downlink multiuser mmWave system with the knowledge of long-term CSI and blockage probability. Specifically, our optimization objective in this work is to minimize the maximum outage probability of all the users. Different from the sum outage probability minimization problem in [START_REF] Iimori | Stochastic learning robust beamforming for millimeter-wave systems with path blockage[END_REF], [START_REF] Zhou | Stochastic learning-based robust beamforming design for RIS-aided millimeter-wave systems in the presence of random blockages[END_REF], the min-max outage probability objective can ensure the quality of service (QoS) performance for the worst-case user in the downlink multiuser system. Due to the non-differentiable objective function, the stochastic gradient descent (SGD) method adopted in [START_REF] Iimori | Stochastic learning robust beamforming for millimeter-wave systems with path blockage[END_REF], [START_REF] Zhou | Stochastic learning-based robust beamforming design for RIS-aided millimeter-wave systems in the presence of random blockages[END_REF] cannot be directly applied. To resolve this problem, two more general and powerful stochastic optimization frameworks are adopted to jointly optimize the active precoding at the BS and the passive beamforming at the RISs.

• The single user case is firstly considered to obtain insights for the robust beamforming design by minimizing the outage probability in the presence of long-term CSI and random blockages. Due to the non-deterministic expression of the probability function in the objective function, we replace it with an expectation function over the long-term CSI and blockage probability. Then, the inner step function in the expectation is approximated by a smooth function which is twice differentiable for the precoding at the BS and the passive beamforming at the RISs, respectively. The resulting expectation optimization problem is solved by adopting the stochastic majorization-minimization (SMM) method in which an upper bound surrogate function of the original differentiable function is constructed for a new channel realization at each iteration. The constructed surrogate problem has closedform solutions and is computationally efficient. We prove that the proposed SMM method is guaranteed to converge to the set of stationary points of the original expectation minimization problem.

• The robust beamforming for a more general multiuser case is designed by solving a maxi-mum outage probability minimization problem. To tackle the non-differentiability of the max objective function, we firstly replace it with its log-sum-exp upper bound. The stochastic successive convex approximation (SSCA) method is then adopted, which has more flexibility in the choice of the surrogate function than the SMM method, and generates closed-form solutions at each iteration. Furthermore, we prove that the final solution generated by the iterative algorithm is guaranteed to converge to the stationary point of the original expectation minimization problem.

• We demonstrate through numerical results that the proposed robust beamforming algorithm outperforms its non-robust counterpart and the robust beamforming in conventional RIS-free systems both in terms of maximum outage probability and minimum effective rate if the blockage probability is high. Moreover, deploying multiple small-size RISs is shown to be more effective than deploying a single large-size RIS in terms of improving the performance of the worst-case user.

The remainder of this paper is organized as follows. Section II introduces the system model.

The outage probability minimization problem is formulated for a single-user case in Section III.

Section IV further investigates the min-max outage probability problem for the multiuser system. Finally, Sections VI and VII report the numerical results and conclusions, respectively.

Notations:

The following mathematical notations and symbols are used throughout this paper.

Vectors and matrices are denoted by boldface lowercase letters and boldface uppercase letters, respectively. The symbols X * , X T , X H , and ||X|| F denote the conjugate, transpose, Hermitian (conjugate transpose), Frobenius norm of matrix X, respectively. The symbol ||x|| 2 denotes 2norm of vector x. The symbols Tr{•}, Re{•}, |•|, λ(•), and (•) denote the trace, real part, modulus, eigenvalue, and angle of a complex number, respectively. diag(x) is a diagonal matrix with the entries of x on its main diagonal. [x] m means the m-th element of the vector x. The Kronecker product between two matrices X and Y is denoted by X ⊗ Y. X Ymeans that X -Y is positive semidefinite. Additionally, the symbol C denotes complex field, R represents real field, and j √ -1 is the imaginary unit. The inner product •, • :

C M ×N × C M ×N → R is defined as X, Y = R{Tr{X H Y}}.

II. SYSTEM MODEL

A. Signal Model

As shown in Fig. 1, we consider an RIS-aided downlink mmWave communication system. In order to ensure high QoS for the users in the presence of random blockages, U RISs, each of which has M passive antennas, are deployed to assist the communication from the BS equipped with N active antennas to K single-antenna users (denoted by K {1, . . . , K}). The RISs are assumed to be connected to controllers that exchange control information with the BS through dedicated channels [START_REF] Pan | Intelligent reflecting surface aided MIMO broadcasting for simultaneous wireless information and power transfer[END_REF], [START_REF] Pan | Multicell MIMO communications relying on intelligent reflecting surfaces[END_REF]. The baseband transmitted signal at the BS is x = Fs, where s ∈ C K×1 ∼ CN(0, I) is the Gaussian data symbol vector and F = [f 1 , . . . , f K ] ∈ C N ×K denotes the full-digital beamforming matrix. The baseband transmit power is limited to the total transmit power P max . Hence, F belongs to the set

S f = {F | ||F|| 2 F ≤ P max }. Let h b,k ∈ C N ×1 , H u ∈ C M ×N and h u,k ∈ C M ×1
denote the channels of the BS-user k, BS-RIS u, and RIS u-user k links, respectively. Then, the received signal intended to the k-th user is expressed as

y k = h H b,k + U u=1 h H u,k E u H u x + n k , (1)
where n k ∼ CN(0, σ 2 k ) is the additive white Gaussian noise (AWGN), and E u = ζdiag([e (u-1)M +1 , . . . , e uM ]) is the reflection coefficient matrix (also known as the passive beamforming matrix) of the u-th RIS. Element e (u-1)M +m is the unit modulus coefficient of the m-th phase shift at the u-th RIS, and ζ ∈ [0, 1] indicates the reflection efficiency. Here, it is assumed that ζ = 1 for investigation of the performance upper bound of RIS. ×N between the BS and the k-th user. The corresponding equivalent reflection coefficient vector is given by e = [e 1 , . . . , e U M +1 ] T ∈ C (U M +1)×1 which belongs to the set S e = {e | |e m | 2 = 1, 1 ≤ m ≤ UM, e U M +1 = 1}. Then, [START_REF] Rappaport | Millimeter wave mobile communications for 5G cellular: it will work![END_REF] can be rewritten into a compact form as

Furthermore, after defining the compact matrices h

k = [h H 1,k , . . . , h H U,k ] H and H = [H H 1 , . . . , H H U ] H , we obtain the equivalent channel G k =   diag(h H k )H h H b,k   ∈ C (U M +1)
y k = e H G k Fs + n k , (2)
and the corresponding achievable signal-to-interference-plus-noise ratio (SINR) is

Γ k (F, e) = |e H G k f k | 2 K i =k |e H G k f i | 2 +σ 2 k . (3)

B. Channel Model

It is important to mention that the perfect instantaneous CSI of a mmWave system is difficult to obtain due to the passive implementation of the RIS reflecting elements. More specifically, estimating the full instantaneous CSI for all links would inevitably require a large training overhead. In this work, since the channel modeling of the RIS-related links is the same as the traditional wireless channel modeling [START_REF] Shen | Modeling and architecture design of intelligent reflecting surfaces using scattering parameter network analysis[END_REF], we adopt a Saleh-Valenzuela (SV) channel model [START_REF] Saleh | A statistical model for indoor multipath propagation[END_REF] to characterize the mmWave channels where only the large-scale fading characteristics is required during the transmission design. In particular, it is assumed that the uniform planar array (UPA) is deployed on the BS and each RIS. The steering vector of the UPA is a (ϕ, φ) in which ϕ(φ) denotes the azimuth (elevation) AoD and the angle-of-arrival (AoA) of the transceivers.

We assume that there are L b,k , L u,k and L b,u sparse scatterers for the BS-user k, the RIS u-user k and the BS-RIS u links, respectively, and each of the scatterers comprises I subpaths. Thus, the mmWave channels can be expressed as

h b,k = g b,k 0 a ϕ AoD b,k,0 , φ AoD b,k,0 + 1 IL b,k L b,k l=1 I i=1 g b,k l,i a ϕ AoD b,k,l,i , φ AoD b,k,l,i , ∀k, (4)
h u,k = g u,k 0 a ϕ AoD u,k,0 , φ AoD u,k,0 + 1 IL u,k L u,k l=1 I i=1 g u,k l,i a ϕ AoD u,k,l,i , φ AoD u,k,l,i , ∀k, ∀u, (5)
H u = g b,u 0 a ϕ AoA u,0 , φ AoA u,0 a ϕ AoD b,0 , φ AoD b,0 H + 1 IL b,u L b,u l=1 I i=1 g b,u l,i a ϕ AoA u,l,i , φ AoA u,l,i a ϕ AoD b,l,i , φ AoD b,l,i H , ∀u, (6)
where, by denoting an arbitrary element q ∈ {(b, k), (u, k), (b, u)} ∀k,∀u , then g q 0 a ϕ AoD q,0 , φ AoD q,0 is the line-of-sight (LoS) component with fading coefficient following the distribution g q 0 ∼ CN(0, ζ q 0 10 PL 10 ), where ζ q 0 = κ 1+κ is the power fraction that corresponds to the Rician factor κ, and PL is the distance-dependent large-scale fading model. The remaining paths are the line-ofsight (NLoS) components whose fading coefficients follow the distribution g q l,i ∼ CN(0, ζ q l 10 PL 10 )

with power fraction ζ q l = 1 (Lq-1)(1+κ) . We assume that the user's locations are quasi-static over milliseconds or even seconds. Therefore, the large-scale fading characteristic parameters, such as the distance-dependent path loss, the numbers of clusters, their power fraction, the cluster central angles and angular beamspreads, change relatively slowly [START_REF] Tse | Fundamentals of Wireless Communication[END_REF] and can be perfectly known by the BS. However, the instantaneous CSI, given by {h b,k , h u,k , H u } 1 , vary during the transmission due to the rapidly varying smallscale fading coefficients {g q 0 , g q l,i }, AoDs and AoAs according to an ergodic stationary process. These AoDs and AoAs can be generated according to a Gaussian distribution, whose expectation values are the cluster central angles and variance values are the angular spreads [START_REF] Tse | Fundamentals of Wireless Communication[END_REF].

Based on the fact that the communication links in the mmWave frequency band are sensitive to the presence of blockages, most existing contributions considered a worst-case scenario where the BS-user links are completely blocked by the obstacles during the whole transmission, while the RIS-related links are not affected by blockages since the locations of the RISs can be appropriately chosen to ensure line-of-sight transmission. However, this assumption may be impractical in many scenarios. Traditionally, blockage effects are incorporated into the shadowing model, along with reflections, scattering, and diffraction [START_REF] Bai | Analysis of blockage effects on urban cellular networks[END_REF]. In contrast, we adopt a recently proposed probabilistic model [START_REF] Iimori | Stochastic learning robust beamforming for millimeter-wave systems with path blockage[END_REF] to characterize the channel uncertainties caused by the presence of random blockages. In particular, by introducing the blockage parameters ω k,l ∈ {0, 1}, 0 ≤ l ≤ L b,k , ∀k ∈ K, which are random variables following a Bernoulli distribution with a blockage probability of p k , the channels in (4) for the BS-user links are modified to

h b,k =ω k,0 g b,k 0 a ϕ AoD b,k,0 , φ AoD b,k,0 + 1 I L b,k l=1 ω k,l I i=1 g b,k l,i a ϕ AoD b,k,l,i , φ AoD b,k,l,i , ∀k, (7)
Analytical results show that the probability of link blocking increases exponentially with the link length [START_REF] Bai | Analysis of blockage effects on urban cellular networks[END_REF], thus we adopt a distance-dependent blockage probability model given by p k (d k ) =

1 With the knowledge of long-term large-scale fading characteristic parameters, each sample of {h b,k , h u,k , Hu} can be obtained by generating {g q 0 , g q l,i }, AoDs and AoAs according to their Gaussian distribution.

max(0, 1 -e -aoutd k +bout ), where physical distance parameter d k between the BS and user k can be extracted from the estimated CSI by the BS, and parameters a out and b out are fit from the data [START_REF] Akdeniz | Millimeter wave channel modeling and cellular capacity evaluation[END_REF].

C. Problem Formulation

To account for the presence of random blockages and to enhance the QoS of the users, we design a robust beamforming by minimizing the maximum outage probability among all users over the random CSI. The formulated optimization problem is given by min

F,e max k∈K Pr{Γ k (F, e) ≤ γ k } (8a) s.t. F ∈ S f (8b) e ∈ S e , (8c)
where the outage probability Pr{Γ k (F, e) ≤ γ k } is the probability that the SINR Γ k (F, e) of the user k is less than the SINR reliability threshold γ k for all possible realizations of the random

channel G = [G 1 , . . . , G K ].
Compared with the sum outage probability minimization problem in [START_REF] Zhou | Stochastic learning-based robust beamforming design for RIS-aided millimeter-wave systems in the presence of random blockages[END_REF], the objective function in [START_REF] Alrabeiah | Deep learning for mmWave beam and blockage prediction using sub-6 GHz channels[END_REF] can ensure the fairness among the users. However, due to the min-max function, the objective function is not smooth or differentiable, and the algorithms in [START_REF] Zhou | Stochastic learning-based robust beamforming design for RIS-aided millimeter-wave systems in the presence of random blockages[END_REF] cannot be directly applied.

III. SINGLE-USER SYSTEM

In this section, we consider a single-user system model in order to obtain some design insights.

By setting K = 1 and dropping the user index, Problem ( 8) is simplified to

min f ,e Pr{Γ (f, e) ≤ γ} (9a) s.t. f ∈ S f (9b) e ∈ S e . (9c)

A. Problem Reformulation

The probability Pr{Γ (f, e) ≤ γ} has no closed-form expression and thus the problem in ( 9) is prohibitively challenging to be solved directly. An alternative solution is to replace the probability function with an equivalent expectation function, i.e., Pr{Γ(f, e)

≤ γ} = E G [I Γ≤γ ]
where I Γ≤γ is the step function of the event Γ ≤ γ. By doing so, various stochastic programming techniques can be used. However, the step function is discontinuous, and the existing stochastic programming methods cannot be directly applied.

To resolve this issue, we approximate the step function with a smooth approximating function

u (x) = 1 1 + e -θx , (10)
where x = γ -Γ and θ is the smooth parameter controlling the approximation error.

By defining f (f, e|G) = u γσ 2 -|e H Gf| 2 , a convenient approximation of Problem ( 9) is

min f ∈S f ,e∈Se g(f, e|G) = E [f (f, e|G)] . (11)

B. Stochastic Majorization-Minimization Method

A simple approach for solving the above problem is the sample average approximation (SAA) method. However, the SAA method is computationally prohibitive since it requires large-sized memory storage due to the fact that the solution obtained at each iteration is calculated by averaging over a large number of channel realizations. To overcome these difficulties, we adopt the widely used SMM [START_REF]Stochastic majorization-minimization algorithms for largescale optimization[END_REF] (also known as stochastic successive minimization [START_REF] Razaviyayn | A stochastic successive minimization method for nonsmooth nonconvex optimization with applications to transceiver design in wireless communication networks[END_REF]) method in which a well-chosen upper bound approximation of the function f (f, e|G) for a new channel sample is constructed at each iteration and the solution is obtained on the average of the new channel sample and the others generated in previous iterations.

The standard of constructing the upper bound approximation function of f (f, e|G) is to make the corresponding surrogate problem easy to solve and even obtain the closed-form solutions. In particular, denote x ∈ {f, e} that belongs to S x ∈ {S f , S e }, the surrogate function f (x,

x i-1 |G)
of f (x|G) around any feasible point x i-1 needs to satisfy the following assumptions [START_REF] Razaviyayn | A stochastic successive minimization method for nonsmooth nonconvex optimization with applications to transceiver design in wireless communication networks[END_REF].

Assumption A (A1) : f x, x i-1 |G is continuous in x for ∀x i-1 ∈ S x . (A2) : f (x i-1 , x i-1 |G) = f (x i-1 |G), ∀x i-1 ∈ S x . (A3) : f (x, x i-1 |G) ≥ f (x|G), ∀x, x i-1 ∈ S x . (A4) : f ′ (x i-1 |G; d) = f ′ (x i-1 |G; d)
, for all x i-1 ∈ S x and feasible directions ∀d at x i-1 .

f ′ (x i-1 |G; d) defines the directional derivative of f (x i-1 |G) in the direction d and is given by

f ′ (x i-1 |G; d) = lim λ→0 f (x i-1 + λd|G) -f (x i-1 |G) λ .
The Assumptions (A2)-(A3) indicate that the surrogate function f (x, x i-1 |G) is a locally upper bound of the original function f (x|G) around the feasible point x i-1 . Assumption (A4) is a derivative consistency condition. To ensure the convergence of the SMM algorithm, we further make the following assumptions [START_REF] Razaviyayn | A stochastic successive minimization method for nonsmooth nonconvex optimization with applications to transceiver design in wireless communication networks[END_REF].

Assumption B

(B1) : The feasible set S x and channel realizations are bounded.

(B2) : The functions f (x, x i-1 |G) and f (x|G), their derivatives, and their second-order derivatives are uniformly bounded.

Since the variables f and e are highly coupled with each other, we adopt an alternating optimization (AO) method to update them. Based on the above assumptions, the variables f and e are updated by solving the following two SMM subproblems:

f n = arg min f ∈S f 1 n n i=1 f f, f i-1 |G i , (12)
and

e n = arg min e∈Se 1 n n i=1 f e, e i-1 |G i . (13)
Here, G 1 , G 2 , ... are some independent samples of the random equivalent channel G. f (f,

f i-1 |G i )
is a surrogate function corresponding to f when e is given, while f (e, e i-1 |G i ) is the corresponding surrogate function of e with given f.

1) Optimizing f: First, we construct f (f, f i-1 |G i ) satisfying Assumptions A and B, which is shown in the following lemma.

Lemma 1 For the twice differentiable function f (f|G i ), we construct its second-order upper bound approximation around any fixed f i-1 , which is given by

f (f, f i-1 |G i ) = 2Re d i,H f f + α i f ||f|| 2 2 +const i f , (14)
where

d i f = m i f -α i f f i-1 , (15a)
m i f = -θe -θx i 1 + e -θx i 2 G i,H e i-1 e i-1,H G i f i-1 , (15b)
x i = γσ 2 -|e i-1,H G i f i-1 | 2 , (15c)
α i f = θ 2 2 P max |e i-1,H G i G i,H e i-1 | 2 , ( 15d
)
const i f = f (f i-1 |G i ) + α i f ||f i-1 || 2 2 -2Re m i,H f f i-1 . (15e)
Proof: See Appendix A.

By using [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF] and ignoring constants, the subproblem in [START_REF] Renzo | Smart radio environments empowered by reconfigurable AI meta-surfaces: An idea whose time has come[END_REF] for updating f is formulated as

min f ∈S f 2Re 1 n n i=1 d i,H f f + 1 n n i=1 α i f ||f|| 2 2 . (16)
Problem ( 16) is convex and can be solved by computing its Lagrange function given by

L(f,κ) = 2Re 1 n n i=1 d i,H f f + 1 n n i=1 α i f ||f|| 2 2 +κ ||f|| 2 2 -P max , (17)
where κ ≥ 0 is a Lagrange multiplier associated with the power constraint. By setting ∂L(f)/∂f * = 0, the globally optimal solution of f at the n-th iteration is derived as

f n = -1 κ + 1 n n i=1 α i f 1 n n i=1 d i f . (18)
Also, (18) must satisfy the power constraint, which yields

|| 1 n n i=1 d i f || 2 2 (κ + 1 n n i=1 α i f ) 2 ≤ P max . (19)
Based on the fact that the left hand side of ( 19) is a decreasing function of κ, we obtain the following closed-form solution

f n =      -1 n i=1 α i f n i=1 d i f , if || n i=1 d i f || 2 2 ( n i=1 α i f ) 2 ≤ P max , - Pmax || n i=1 d i f || 2 2 n i=1 d i f , otherwise. (20)
The first option in ( 20) is based on κ = 0. The second option is due to the fact that there must exist a κ > 0 that (19) holds with equality.

2) Optimizing e: As for the update of e with given f, we first construct a surrogate function corresponding to e in the following lemma.

Lemma 2 For the twice differentiable function f (e|G i ), we construct its second-order upper bound approximation around any feasible e i-1 , which is given by

f (e, e i-1 |G i ) = 2Re d i,H e e + const i e , (21)
where

d i e = m i e -α i e e i-1 , (22a)
m i e = -θe -θx i 1 + e -θx i 2 G i f i-1 f i-1,H G i,H e i-1 , (22b)
α i e = θ 2 2 (UM + 1)|f i-1,H G i,H G i f i-1 | 2 , ( 22c
)
const i e = f (e i-1 |G i ) + 2(UM + 1)α i e -2Re m i,H e e i-1 . ( 22d
)
Proof: The proof of Lemma 2 is similar to that of Lemma 1 and hence omitted for brevity.

By substituting [START_REF] Zhi | Power scaling law analysis and phase shift optimization of RIS-aided massive MIMO systems with statistical CSI[END_REF] into the objective function of subproblem [START_REF]Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and the road ahead[END_REF] and ignoring constants, we obtain:

min e∈Se 2Re 1 n n i=1 d i,H e e . ( 23
)
The globally optimal solution of the above problem is given by

e n = exp j n i=1 d i e / n i=1 d i e U M +1 , (24)
where [•] m means the m-th element of vector, j √ -1 is the imaginary unit, (•) denotes the angle of a complex number, and exp {j (•)} is an element-wise operation.

C. Algorithm Development

Under the SMM framework, closed-form solutions of f in [START_REF] Zhou | Robust beamforming design for intelligent reflecting surface aided MISO communication systems[END_REF] and of e in (24) at each iteration are obtained. Such simple closed-form solutions can greatly reduce the computational complexity. Algorithm 1 summarizes the proposed robust beamforming design based on the SMM-based outage probability minimum problem for RIS-aided single-user mmWave systems in which the BS-user links experience random blockages. The proposed algorithm is referred to as SMM-OutMin. Note that the convergence speed of the SMM-based algorithm might be affected by the tightness of the upper bounds in Lemma 1 and Lemma 2, thus SQUAREM [START_REF] Varadhan | Simple and globally convergent methods for accelerating the convergence of any EM algorithm[END_REF] is adopted to accelerate the SMM-based algorithm. Denote by F (f n ) and F (e n ) the objective function values of Problem [START_REF] Pan | Reconfigurable intelligent surface for 6G and beyond: Motivations, principles, applications, and research directions[END_REF] and Problem [START_REF] Wu | Joint active and passive beamforming optimization for intelligent reflecting surface assisted SWIPT under QoS constraints[END_REF] in the n-th iteration, respectively. P S (•) is an operation mapped to nonlinear constraint set. For the power constraint set S f , P S (•) can be

P S (x) = (x) ||x|| 2 ||f 2 || 2 .
For the unit-modulus constraint set S e , P S (•) can be element-wise operation, i.e., P S (x) = j (x). Steps 10 to 13 and steps 21 to 24 are used to maintain the convergence property of the objective function values.

Algorithm 1 SMM-OutMin Algorithm

Initialize: Initialize f 0 ∈ S f and e 0 ∈ S e . Set n = 1.

1: repeat 2:
Obtain the sample channel G n .

3:

Set e = e n-1 .

4:

Obtain f 1 according to [START_REF] Zhou | Robust beamforming design for intelligent reflecting surface aided MISO communication systems[END_REF] based on f n-1 .

5:

Obtain f 2 according to [START_REF] Zhou | Robust beamforming design for intelligent reflecting surface aided MISO communication systems[END_REF] based on f 1 . 6:

j 1 = f 1 -f n-1 . 7: j 2 = f 2 -f 1 -j 1 . 8: ω = -||j 1 || 2 ||j 2 || 2 .
9:

f n = -P F (f n-1 -2ωj 1 + ω 2 j 2 ). 10: while |F (f n ) -F (f n-1 )|≤ |F (f n-1 ) -F (f n-2 )| do 11: ω = (ω -1)/2.
12:

f n = -P S (f n-1 -2ωj 1 + ω 2 j 2 ).
13:

end while 14:

Set f = f n .

15:

Obtain e 1 according to [START_REF] Hong | Artificial-noise-aided secure mimo wireless communications via intelligent reflecting surface[END_REF] based on e n-1 .

16:

Obtain e 2 according to [START_REF] Hong | Artificial-noise-aided secure mimo wireless communications via intelligent reflecting surface[END_REF] based on e 1 . 17:

j 1 = e 1 -e n-1 .
18:

j 2 = e 2 -e 1 -j 1 .
19:

ω = -||j 1 || 2 ||j 2 || 2 .
20:

e n = -P S (e n-1 -2ωj 1 + ω 2 j 2 ). ω = (ω -1)/2.

23:

e n = -P S (e n-1 -2ωj 1 + ω 2 j 2 ).

24:

end while 25:

n = n + 1. 26: until ||F (f n ) -F (f n-1 ) || 2 → 0 and ||F (e n ) -F (e n-1 ) || 2 → 0.
1) Convergence analysis: The convergence of Algorithm 1 is given in the following theorem.

Theorem 1 Suppose Assumptions A and B are satisfied. Then the sequence of the solutions obtained in each iteration of Algorithm 1 converge to the set of stationary points of Problem [START_REF] Zhou | Stochastic learning-based robust beamforming design for RIS-aided millimeter-wave systems in the presence of random blockages[END_REF] almost surely.

Proof: See Appendix B.

2) Complexity analysis: The computational complexity for updating f n and e n at each iteration mainly depends on the computation of ( 20) and ( 24), respectively. In particular, due to the update rule in

{ n i=1 α i f , n i=1 d i f , n i=1 d i e }, only {α n e , d n f ,
d n e } needs to be calculated at the n-th iteration. Therefore, the approximate complexity of each iteration is given by O(4UMN +12N).

IV. MULTIUSER SYSTEM

In this section, we consider the general multiuser setup and solve Problem [START_REF] Alrabeiah | Deep learning for mmWave beam and blockage prediction using sub-6 GHz channels[END_REF]. Problem ( 8) is more challenging than problem (9) due to the non-differentiable objective function. Furthermore, the complex objective function complicates the use of the SMM algorithm, thus we extend the SMM method to a general algorithm to solve Problem [START_REF] Alrabeiah | Deep learning for mmWave beam and blockage prediction using sub-6 GHz channels[END_REF].

A. Problem Reformulation

We first approximate the probability function in the original Problem (8) by using the expectation of the smooth function in [START_REF] Iimori | Stochastic learning robust beamforming for millimeter-wave systems with path blockage[END_REF]. By defining

f k (F, e|G) = u e H G k FΥ k F H G H k e + γ k σ 2 k , in which Υ k is a diagonal matrix whose diagonal entry is γ k except in the k-th diagonal element that is -1, the approximate objective function is max k∈K E [f k (F, e|G)].
However, the obtained function is still intractable since the maximization operation couples f k , ∀k and the different channel states due to the expectation operation. This issue motivates us to use the following Jensen inequality

max k∈K E [f k (F, e|G)] ≤ E max k∈K f k (F, e|G) , (25)
due to the fact that the max function max k∈K {x 1 , . . . , x K } is convex [START_REF] Boyd | Convex optimization[END_REF].

Furthermore, the non-differentiable max function, max k∈K f k (F, e|G), is approximated by adopting a smooth log-sum-exp upper-bound [START_REF] Xu | Smoothing method for minimax problems[END_REF] max

k∈K f k (F, e|G) ≈ F (F, e|G) = µ ln k∈K exp 1 µ f k (F, e|G) , (26)
where µ > 0 is a smoothing parameter satisfying

max k∈K f k (F, e|G) ≤ F (F, e|G) ≤ max k∈K f k (F, e|G) + 1 µ log (|K|) . ( 27
)
When µ is chosen appropriately, the smooth approximation of Problem ( 8) is approximately reformulated as

min F∈S f ,e∈Se G (F, e|G) = E [F (F, e|G)] . (28)

B. Stochastic Successive Convex Approximation Method

Similar to Problem [START_REF] Zhou | Stochastic learning-based robust beamforming design for RIS-aided millimeter-wave systems in the presence of random blockages[END_REF], Problem [START_REF] Perović | Channel capacity optimization using reconfigurable intelligent surfaces in indoor mmWave environments[END_REF] can still be solved by adopting the SMM method.

However, the function F (F, e|G) in ( 26) is much more complex and its second-order derivative, which is necessary to construct the upper bound surrogate function of F (F, e|G) as shown in Appendix A, is not easy to be calculated. Furthermore, the coefficient of the second-order term in the final upper-bound surrogate function of F (F, e|G) (α i f in ( 14)) can be quite loose, resulting in a very slow convergence rate of the SMM algorithm. Therefore, in this section, we adopt a flexible SSCA method to address the above issues. The surrogate functions employed by the SSCA method do not need to be an upper bound of the original function but they need only to preserve the first-order property of the original function.

Accordingly, the surrogate functions need to satisfy Assumption B and the following assumptions [START_REF] Razaviyayn | Successive convex approximation: Analysis and applications[END_REF].

Assumption C (C1) : F x, x i-1 |G is strongly convex in x for ∀x i-1 ∈ S x . (C2) : F (x i-1 , x i-1 |G) = F (x i-1 |G), ∀x i-1 ∈ S x . (C3) : ∇ x F (x i-1 , x i-1 |G) = ∇ x F (x i-1 |G), ∀x, x i-1 ∈ S x .
Assumption C cannot ensure that the sequences of the approximate objective values are monotonically decreasing at each iteration. Nevertheless, to guarantee convergence, the variables can be updated by choosing an appropriate step size at each iteration that yields a sufficient decrease of the objective value. Based on the above assumptions, we choose the proximal gradient-like approximation to construct the surrogate function, which is

F x, x i-1 |G =F x i-1 |G + ∇ x F (x i-1 |G) T (x -x i-1 ) + τ i 2 ||x -x i-1 || 2 , ( 29
)
where τ i can be any positive number.

1) Optimizing f: By using ( 29), we construct a surrogate function as

F F, F i-1 |G = F F i-1 |G + Tr ∇ F F F i-1 |G T (F -F i-1 ) + Tr ∇ F * F F i-1 |G T (F * -F i-1, * ) + τ i 2 ||F -F i-1 || 2 F = F F i-1 |G + 2 k∈K l i k Re Tr Υ k F i-1,H G H k e i-1 e i-1,H G k (F -F i-1 ) + τ i 2 ||F -F i-1 || 2 F = 2Re Tr P i,H f F + τ i 2 ||F|| 2 F +cons1 i , ( 30
)
for F around the fixed F i-1 when e is given.

The parameters in [START_REF] Ning | Channel estimation and transmission for intelligent reflecting surface assisted THz communications[END_REF] are as follows

P i f = W i f - τ i 2 F i-1 , (31a)
W i f = k∈K l i k G i,H k e i-1 e i-1,H G i k F i-1 Υ k , ( 31b
)
l i k = exp 1 µ f k (F i-1 , e i-1 |G i ) k∈K exp 1 µ f k (F i-1 , e i-1 |G i ) θe -θx i k 1 + e -θx i k 2 , (31c)
x i k = e i-1,H G i k F i-1 Υ k F i-1,H G i,H k e i-1 + γ k σ 2 k , ( 31d
)
cons1 i = F F i-1 |G i + τ i 2 ||F i-1 || 2 F -2Re Tr W i,H f F i-1 . ( 31e
)
By using [START_REF] Ning | Channel estimation and transmission for intelligent reflecting surface assisted THz communications[END_REF], the subproblem of Problem [START_REF] Perović | Channel capacity optimization using reconfigurable intelligent surfaces in indoor mmWave environments[END_REF] corresponding to F at the n-th iteration is formulated as

min F∈S F 1 n n i=1 F F, F i-1 |G i . ( 32
)
The method to solve Problem [START_REF] Shen | Modeling and architecture design of intelligent reflecting surfaces using scattering parameter network analysis[END_REF] is the same as that of solving Problem ( 16), thus we directly provide the global minimizer of Problem [START_REF] Shen | Modeling and architecture design of intelligent reflecting surfaces using scattering parameter network analysis[END_REF] as

F n =      -2 n i=1 τ i n i=1 P i f , if 4|| n i=1 P i f || 2 F ( n i=1 τ i ) 2 ≤ P max , - Pmax || n i=1 P i f || 2 F n i=1 P i f , otherwise. ( 33
)
2) Optimizing e: Furthermore, with fixed F, the subproblem of Problem [START_REF] Perović | Channel capacity optimization using reconfigurable intelligent surfaces in indoor mmWave environments[END_REF] corresponding to e at the n-th iteration is

min e∈Se 1 n n i=1 F e, e i-1 |G i , (34)
where F (e, e i-1 |G i ) = 2Re p i,H e e + cons2 i , and

p i e = w i e - τ i 2 e i-1 , (35a)
w i e = k∈K l k G k F i-1 Υ k F i-1,H G H e i-1 , ( 35b
) cons2 i = F e i-1 |G i + τ i (UM + 1) -2Re w i,H e e i-1 . (35c)
Finally, the minimizer of Problem ( 34) is

e n = exp j n i=1 p i e / n i=1 p i e U M +1 . ( 36
)

C. Algorithm development

The simple closed-form solutions of f in [START_REF] Saleh | A statistical model for indoor multipath propagation[END_REF] and of e in [START_REF] Razaviyayn | A stochastic successive minimization method for nonsmooth nonconvex optimization with applications to transceiver design in wireless communication networks[END_REF] can greatly reduce the computational complexity. Algorithm 2 summarizes the proposed SSCA-based robust beamforming design for RIS-aided multiuser mmWave systems in which the BS-user links experience random blockages. The proposed algorithm is referred to as SSCA-OutMin.

Algorithm 2 SSCA-OutMin Algorithm

Initialize: Initialize F 0 ∈ S f and e 0 ∈ S e . Set n = 0.

1: repeat 2:

n = n + 1.

3:

Obtain the sample channel G n .

4:

Calculate F n according to [START_REF] Saleh | A statistical model for indoor multipath propagation[END_REF].

5:

Update

F n = F n-1 + ξ n f F n -F n-1 . 6:
Calculate e n according to [START_REF] Razaviyayn | A stochastic successive minimization method for nonsmooth nonconvex optimization with applications to transceiver design in wireless communication networks[END_REF].

7:

Update e n = e n-1 + ξ n e ( e ne n-1 ).

8: until ||F n -F n-1 || 2 F → 0 and ||e n -e n-1 || 2 → 0.

1)

Step-size selection: It is worth noting that the approximation in [START_REF] Ning | Channel estimation and transmission for intelligent reflecting surface assisted THz communications[END_REF] has the same form as that in [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF]. However, τ i in the SSCA method can be any positive number, and F (F, F i-1 |G) might no longer be a global upper bound of F (F|G). In this case, the step sizes ξ n f and ξ n e need to be carefully chosen to ensure convergence.

As an example, we take ξ n f to illustrate the update rule, which is a line-search (also called Armijo step-size) rule: Consider ξ 0 f > 0 and c 1,f , c 2,f ∈ (0, 1). Let ξ n f be the largest element in {ξ 0 f c t 2,f } t=0,1,... such that

F F n-1 +ξ n f F n -F n-1 ≤ F (F n-1 ) + c 1,f ξ n f Tr ∇ F F (F n-1 ) T F n -F n-1 . ( 37
)
Theorem 2 If {ξ n f } n=1,2,... is chosen according to the line-search rule, then

lim n→∞ || F n -F n-1 || = 0.
Proof: See Theorem 7 in [START_REF] Dunford | Successive convex approximation: Analysis and applications[END_REF].

2) Convergence analysis: The convergence of Algorithm 2 is given in the following theorem.

Theorem 3 Suppose Assumptions B and C are satisfied. Then every limit point of the iterations generated by Algorithm 2 is a stationary point of Problem ( 28) almost surely.

Proof: See Appendix C.

3) Complexity analysis: The computational complexity for updating f n and e n at each iteration mainly depend on the computation of ( 33) and [START_REF] Razaviyayn | A stochastic successive minimization method for nonsmooth nonconvex optimization with applications to transceiver design in wireless communication networks[END_REF], respectively. In particular, only {P n f , p n e } needs to be calculated at the n-th iteration. Therefore, the approximate complexity of each iteration is given by O((K + 2)2UMN + UMK + NK + (N + 2)K 2 + 2N).

H G 0 k || 2 2 . (38)
Problem ( 38) can be efficiently solved by using the SDR method as follows

max E t (39a) s.t. Tr{G 0 k G 0,H k E} ≥ t, ∀k ∈ K (39b) E 0, rank(E) = 1, [E] m,m = 1, ∀m, (39c)
where E = ee H and t is an auxiliary variable.

Furthermore, F is initialized by using the maximum-ratio transmission (MRT) method as

F 0 = P max G 0 e 0 ||G 0 e 0 || . (40)

V. NUMERICAL RESULTS AND DISCUSSION

A. Simulation setup

In this section, we numerically evaluate the performance of the proposed algorithms. All experiments are performed on a PC with a 1.99 GHz i7-8550U CPU and 16 GB RAM. We adopt a polar coordinate to describe the simulated system setup as shown in Fig. 2, where the BS is located at (0 m, 0), and the two RISs are deployed in the locations (50 m, 0) and ( 50 where σ 2 ¸denotes the log-normal shadowing variance. The mmWave system operates on carrier frequency f c = 28 GHz and bandwidth 20 MHz. Since the macro-scattering environment between the BS and the users are complex, only NLoS clusters are assumed to exist in the BS-user links, i.e., the Rician factor is κ = 0. The parameter settings in PL of NLoS are α = 3.5 and σ ξ = 8.2 dB [START_REF]Technical Specification Group Radio Access Network; Study on Channel Model for Frequencies from 0.5 to 100 GHz (Release 16)[END_REF]. In practice, the RISs can be installed such that the BS-RIS links and the RIS-user links are blockage-free. Thus, the channels in ( 5) and ( 6) contain only LoS cluster with a Rician factor κ → ∞. The parameters in PL of LoS are α = 2 and σ ξ = 4 dB according to [START_REF]Technical Specification Group Radio Access Network; Study on Channel Model for Frequencies from 0.5 to 100 GHz (Release 16)[END_REF].

Unless stated otherwise, we assume L b,k = L u,k = L b,u = 5 and I = 20. The transmit power limit of the BS is P max = 30 dBm and the noise power at each user is σ 2 1 = . . . = σ 2 K = -94 dBm. For simplicity, we consider an equal blockage probability, p k,l = p block , ∀k, l, and an equal target SINR, γ = γ 1 = . . . = γ K , which yields the target rate R targ = log 2 (1 + γ). The smooth parameters are chosen to be θ =

1 max ∀k∈K |x 0 k | and µ = 1
100K . To evaluate the performance of the proposed stochastic optimization algorithms, we consider the following benchmark schemes. 1) SMRT: In this scheme, the active precoding is updated by using stochastic maximum-ratio transmission (SMRT), that is

F n+1 = 1 n n i=1 G i,H e i || 1 n n i=1 G i,H e i || 2 √ P max .
The passive beamforming is still updated by using the SMM or SSCA methods. 2) NoRIS: In this case, no RIS is employed and the optimal active precoding is obtained by using the SMM or SSCA methods. 3) No-robust: In this scheme, the beamforming is designed by using the SMM or SSCA methods by taking into account the random small-scale parameters while assuming p block = 0. 4) Imperfect CSI: In this scheme, the beamforming is designed by using the SMM or SSCA methods based on the imperfect CSI of cluster central angles. In specific, we assume that the modulus of the cluster central angle estimation error is 0.01. 5) SAA: In this scheme, we generate 300 independent channel realizations in advance, the solutions at each iteration is the average over these 300 channel samples, and the surrogate function used at each iteration is obtained by adopting the MM or SCA methods. We take the beamforming design in the singleuser case as an example. By modifying Problems ( 12) and ( 13), the beamforming designed by using the SAA-MM method is updated as follows

f n = arg min f ∈S f 1 300 300 i=1 f f, f n-1 |G i , (41)
In order to demonstrate the robustness of the proposed algorithms, we consider two performance metrics: the outage probability and the effective rate. In particular, the outage probability of each user is calculated by averaging over 1000 independent channel realizations. The corresponding effective rate of the k-th user is defined as

R eff,k E[log 2 (1+Γ k (F, e))] if Γ k (F, e) ≥ γ
and R eff 0 otherwise.

B. Convergence

Fig. 3 investigates the convergence behavior of the considered stochastic optimization algorithms. For comparison, we consider a single-user case containing RIS 1 in Fig. 2, and the other parameters are given in Fig. 3. In Fig. 3, the coordinate value on the y-axis is the objective value of Problem ( 16) or [START_REF] Shen | Modeling and architecture design of intelligent reflecting surfaces using scattering parameter network analysis[END_REF], and not the actual outage probability of the original problem. It is observed from Fig. 3 that the SMM and SSCA algorithms are characterized by an oscillatory convergence which depends on the random channel generations at each iteration. On the other hand, using 300 channel realizations for each iteration leads to the monotonic convergence of the SAA algorithm when adopting a monotonically decreasing surrogate function for each channel realization. Although the SAA algorithm requires the least number of iterations to converge, it is much more computationally demanding than the other two algorithms. This fact can be observed in Table I which compares the CPU time consumption of each iteration and total CPU time consumption of iteration convergence for the three considered algorithms. Theoretically, the computational complexity of each iteration of the SAA algorithm is 300 times higher than that of the SMM or SAA algorithms, because at each iteration of the SAA algorithm, each parameter needs to be calculated 300 times for all channel realizations.

C. Single-user system

We consider a single-user system containing RIS 1 in Fig. 2 with the target rate of R targ = 0.5 bps/Hz. Fig. 4 illustrates the performance of different algorithms as a function of the blockage probability. First, it can be seen that SMM-based beamforming in the RIS-aided mmWave system of M = 64 outperforms the NoRIS scheme when p block ≥ 0.1. The main reason is that the direct BS-user channel may be much stronger than the cascaded BS-RIS-user channel as the latter experiences the double path loss effect, which is significant if the path-loss is large. Therefore, when the blockage probability is small (p block ≤ 0.1), the BS tends to allocate the transmit power to the stronger direct path, thus reducing the contribution of the RIS to the system performance. However, by increasing the number of reflecting elements at the RIS to M = 128, our proposed SMM algorithm outperforms the NoRIS system in the whole blockage probability region (i.e., 0 ≤ p block ≤ 1). The reason is that the RIS-aided reflecting channels starts to compensate for the performance loss caused by the blockages when p block = 0. In addition, it can be seen that the performance of the SMM-based beamforming is the same as that of the SAA-based beamforming. Last but not least, compared with the No-robust RIS-aided scenario and the scheme with imperfect CSI, the robust RIS beamforming can significantly compensate for the performance loss caused by the presence of random blockages.

D. Multiuser system

In this section, we consider a multiuser system with K = 3 users, in which the target rate is R targ = 0.1 bps/Hz. As we can see from Fig. 6(a), compared with the No-robust RIS-aided scenario, the robust beamforming in the RIS-aided scheme in a multiuser system can efficiently improve the maximum outage probability. In addition, different from the single-user system of M = 128, the RIS-aided setup outperforms the NoRIS scheme only when p block ≥ 0.2. This fact is consistent with the observation of M = 64 in Fig. 4 (a), that is, the performance gain of an RIS starts to manifest when p block = 0.2. Further increasing the number of RISs and the size of each RIS can reduce the maximum outage probability in the whole blockage probability Specifically, the minimum effective rate of the robust designs is better than that of the non-robust schemes. Moreover, the contribution of the RIS in improving the minimum effective rate is also apparent for the almost blockage probability region (i.e., p block ≤ 0.8). Finally, we observe that multiple RISs with a relatively small size can improve the worst-case user's performance, and outperform the performance offered by a single RIS with large size, especially in the multi-user scenario. This is because the optimal deployment of multiple RISs can provide more spatial diversity for users who are far away from the RIS, thus improving the system performance.

Finally, Fig. 7 investigates the performance as a function of the number of users. For a fair comparison, we consider the setup N = 32, M = 128 and U = 2. Compared with the NoRIS scheme, the contribution of multiple RISs in reducing the maximum outage probability shown in Fig. 7 decreases as the increase of K. When K ≥ 5, the RIS scheme cannot guarantee the QoS performance of the worst-case user due to the high outage probability. Closed-form solutions have been derived at each iteration by adopting the SMM and SSCA

The corresponding directional derivative of f (f) is

-θe -θx i 1 + e -θx i 2 2Re f i-1,H G i,H e i-1 e i-1,H G f -f i-1 , (45)
where x i is given in (15c).

Assumption (A4) is satisfied only when ( 44) and ( 45) are equal, yielding

m i f = - θe -θx i 1 + e -θx i 2 G i,H e i-1 e i-1,H G i f i-1 . (46)
In order for Assumption (A3) to hold, it is sufficient to show that f (f, f i-1 ) is an upper bound for each linear cut in any direction. In particular, let f = f i-1 + ξ( ff i-1 ), ∀ξ ∈ [0, 1], we need to show

f f i-1 + ξ( f -f i-1 ) ≤ f f i-1 + 2ξRe m i,H f ( f -f i-1 ) + ξ 2 ( f -f i-1 ) H M i f ( f -f i-1 ). ( 47
) Define L(ξ) = f (f i-1 + ξ( f -f i-1 )) and l(ξ) = γσ 2 -|e i-1,H G i (f i-1 + ξ( f -f i-1 ))| 2 . ( 47
)
holds if the second-order derivative of L(ξ) is no more than that of the right hand side of (47).

The corresponding sufficient condition can be formulated as

∂ 2 L(ξ) ∂ξ 2 ≤ 2Tr ( f -f i-1 ) H M i f ( f -f i-1 ) . ( 48
)
Before deriving the expression of ∂ 2 L(ξ)/∂ξ 2 , we calculate the first-order derivative of L(ξ), as follows

∂L(ξ) ∂ξ = g(ξ)∇ ξ l(ξ), (49)
where g(ξ) = θe -θl(ξ) (1+e -θl(ξ) ) 2 , ∇ ξ l(ξ) = -2Re{q H ( ff i-1 )}, and q = G i,H e i-1 e i-1,H G i (f i-1 + ξ( ff i-1 )).

Then, the second-order derivative is derived as

∂ 2 L(ξ)
∂ξ 2 = g(ξ)∇ 2 ξ l(ξ) -θg(ξ)∇ ξ l(ξ) (∇ ξ l(ξ)) T + 2 1 + e -θl(ξ) g(ξ)∇ ξ l(ξ) (g(ξ)∇ ξ l(ξ)) T ,

where ∇ 2 ξ l(ξ) = -2Re{( ff i-1 ) H Θ( ff i-1 )} and Θ = ξG i,H e i-1 e i-1,H G i . (50) is rewritten as a quadratic form of t = ff i-1 , as follows For convenience, we select M i f = α i f I = λ max (Φ) I. Then, f (f, f i-1 ) in ( 43) is designed to be

∂ 2 L(ξ) ∂ξ 2 =   t t *   H Φ   t t *   , (51)
f f, f i-1 = f f i-1 + 2Re m i,H f (f -f i-1 ) + α i f ||f -f i-1 || 2 2 = 2Re d i,H f f + α i f ||f|| 2 2 +const i f ,
where d i,H f , α i f and const i f are defined in Lemma 1. The deterministic expression of λ max (Φ) is difficult to obtain, therefore we derive the upper bound shown as follows compact set S x due to the assumption (B2) and the use of mean value theorem. Thus, by restricting to a subsequence, we have (56). Furthermore, the families of functions {ĝ n (x)} are also equicontinuous and bounded over a compact set S x due to the assumption (B2) that ||∇ x f (x, x i-1 , G)|| is bounded. Hence the Arzel à -Ascoli theorem [START_REF] Dunford | Linear operators. Part 1: General theory[END_REF] implies that, by restricting to a subsequence, there exists a uniformly continuous function ĝ(x) such that (57) and (58) hold.

On the other hand, the update rule of Algorithm 1 leads the following lemma.

Lemma 4 lim n→∞ |ĝ n (x n ) -g n (x n )|= 0, almost surely.

Proof: The proof of Lemma 4 is the same with that of ( [START_REF] Razaviyayn | A stochastic successive minimization method for nonsmooth nonconvex optimization with applications to transceiver design in wireless communication networks[END_REF], Lemma 1) and is omitted for conciseness.

Assumption (A3) implies that ĝn j (x) ≥ g n j (x), ∀x ∈ S x . combining it with (55) and (57), we obtain ĝ(x) ≥ g(x), ∀x ∈ S x .

Moreover, combining Lemma 4 with ( 56) and (58), it yields

ĝ(x) = g(x). (60)
Then, (59) and (60) imply that x is a minimizer of function ĝ(x) -g(x), hence the first-order optimality condition satisfies ∇ĝ(x) -∇g(x) = 0.

(61)

Due to the fact that x is the limit point of Problem [START_REF] Renzo | Smart radio environments empowered by reconfigurable AI meta-surfaces: An idea whose time has come[END_REF] or Problem (13), we have ĝ(x) ≤ ĝ(x), ∀x ∈ S x , which implies that ∇ĝ(x), x -x ≥ 0, , ∀x ∈ S x .

Combining this with (61), we obtain

∇g(x), x -x ≥ 0, , ∀x ∈ S x , (63)
which means that the directional derivative of the objective function g(x) is non-negative for every feasible direction at x. Recalling that x ∈ {f, e} and defining the limit points { f, ē}, (63) is equivalent to      ∇g( f), f -f ≥ 0, ∀f ∈ S f , ∇g(ē), e -ē ≥ 0, ∀e ∈ S e .

Since that x ∈ {F, e}, we define the limit points { F, ē}, (72) is then equivalent to      ∇G( F), F -F ≥ 0, ∀F ∈ S f , ∇G(ē), e -ē ≥ 0, ∀e ∈ S e .

Therefore, according to [START_REF] Dunford | Successive convex approximation: Analysis and applications[END_REF], { F, ē} is a stationary point of Problem [START_REF] Perović | Channel capacity optimization using reconfigurable intelligent surfaces in indoor mmWave environments[END_REF] due to the regularity of G(•).

Fig. 1 :

 1 Fig. 1: Multiple RIS-aided mmWave communication system.

21 :

 21 while |F (e n ) -F (e n-1 ) |≤ |F (e n-1 ) -F (e n-2 ) | do 22:

4 )

 4 Initial point: Problem (28) has, in general, multiple local minima points due to the nonconvex unit-modulus constraint and e ∈ S e . The accurate selection of the initial points in Algorithm 2 plays an important role for the convergence speed and the quality of the obtained local solution. To that end, we first initialize e to maximize the minimum equivalent total channel gain, resulting in the following optimization problem e 0 = arg max e∈Se min k∈K ||e

Fig. 2 :

 2 Fig. 2: The simulated system setup.

m, π/ 6 )

 6 which are close to the users. The users are randomly placed in a range with polar diameter d k ∈[50 m, 80 m] and polar angle ϑ ∈ [0, π/6], where d k is used to calculate the distancedependent blockage probability. The large-scale fading, reported in urban micro (UMi)-street canyon scenario [42], is modeled as PL = 32.4 + 20 log 10 (f c ) + 10α log 10 (D) + ξ in dB with link distance D (in meters), path loss exponent α, and log-normal shadowing ξ ∼ CN(0, σ 2 ¸)

  e n-1 |G i .

Fig. 3 :

 3 Fig. 3: Convergence behavior of different algorithms, when N = 8, M = 128, K = 1, U = 1, and R targ = 0.5 bps/Hz.

Fig. 5

 5 Fig.5shows the impact of the size of the RIS and the size of the antenna array at the BS on the outage probability by using the SMM algorithm. It can be observed from Fig.5(a) that when BS is equipped with N = 8 antennas, the RIS plays a significant role in guaranteeing the desired user's QoS and improving the system robustness with increased number of reflection elements (M : 64 → 256). In other words, a large-size RIS with M ≥ 224 can still reduce the outage probability to about 0.1, even if the direct channel is completely blocked. Finally in Fig.5(b) we observe that increasing the number of transmit antennas (N : 8 → 64) has the same ability to reduce the outage probability as increasing the number of reflection elements (M : 64 → 256), except when the direct channel is completely blocked.

Fig. 4 :

 4 Fig. 4: Comparison of outage probability and effective rate as a function of the blockage probability p block for N = 8, K = 1, U = 1, and R targ = 0.5 bps/Hz.

Fig. 5 :

 5 Fig. 5: Comparison of outage probability as a function of M with fixed N = 8 and of N with fixed M = 128, when K = 1, U = 1, and R targ = 0.5 bps/Hz.

Fig. 6 :

 6 Fig. 6: Comparison of maximum outage probability and minimum effective rate as a function of the blockage probability p block for N = 16, K = 3, and R targ = 0.1 bps/Hz.

whereΦ

  = g(ξ) 2 1 + e -θl(ξ) g(ξ) -θ ξ)I 2 ⊗ Θ.Furthermore, we also manipulate the right hand side of (48) into the same form as in (51)

λ 2 2

 2 max (Φ) (p1) ≤ 2 1 + e -θl(ξ) g 2 (ξ)λ max ξ)λ min (I 2 ⊗ Θ) -θg(ξ)λ min = 4 1 + e -θl(ξ) g 2 (ξ)||q||

TABLE I :

 I Comparison of the CPU time

	Algorithms The CPU time (sec) per iteration	The CPU time (sec)
	SMM	0.0025	1.8750
	SSCA	0.0042	4.6719
	SAA	0.3557	20.9844

In this subsection, G i is dropped for simplicity, i.e., f

is twice differentiable, we propose a second-order approximation to upper bound f (f) at any fixed point f i-1 :

where m i f and M i f are to be designed to satisfy Assumption A. Assumptions (A1) and (A2) are readily satisfied. Assumption (A4) is a derivative consistency

The above inequalities are due to the following mathematical properties:

(p1): A and B are Hermitian matrices:

, where the equality holds when l(ξ) = 0.

(p4): A is positive semidefinite with maximum eigenvalue λ max (A) and B is positive semidefinite:

Hence, the proof is completed.

APPENDIX B THE PROOF OF THEOREM 1

Define the random functions

To state the convergence result, we need the following lemmas.

Lemma 3 Suppose Assumptions A and B are satisfied and define a limit point x of a subsequence {x n j } ∞ j=1 , then there exists uniformly continuous functions g(x) and ĝ(x) such that

Proof: First, it follows that f (x, G) is bounded for ∀x ∈ S x and all channel realizations due to the assumption (B2), and therefore, (55) holds by using the strong law of large numbers [START_REF] Fristedt | A Modern Approach to Probability Theory[END_REF]. Also, the families of functions {g n j (x)} are equicontinuous and bounded over a Therefore, according to [START_REF] Dunford | Successive convex approximation: Analysis and applications[END_REF], { f, ē} is a stationary point of Problem [START_REF] Zhou | Stochastic learning-based robust beamforming design for RIS-aided millimeter-wave systems in the presence of random blockages[END_REF] due to the regularity of g(•).

APPENDIX C THE PROOF OF THEOREM 1

Define the random functions

To state the convergence result, we need the following lemmas.

Lemma 5 Suppose Assumptions B and C are satisfied and define a limit point x of a subsequence {x n j } ∞ j=1 , then there exists uniformly continuous functions G(x) and Ĝ(x) such that

Ĝ(x) = lim j→∞ Ĝn j (x n j ).

Proof: The proof of Lemma 5 is the same with that of Lemma 3 and is omitted for brevity.

On the other hand, x n j is the minimizer of Ĝn j (x), thus Ĝn j (x n j ) ≤ Ĝn j (x), ∀x ∈ S x .

(70)

Assuming j → ∞, and combining (68) and (69), we obtain Ĝ(x) ≤ Ĝ(x), ∀x ∈ S x , which implies that its first-order optimality condition satisfies ∇ Ĝ(x), x -x ≥ 0, , ∀x ∈ S x .

By combining (71) and Assumption (C3), we finally obtain ∇G(x), x -x ≥ 0, , ∀x ∈ S x .

(72)