Joint Active-Passive Beamforming and User Association in IRS-Assisted mmWave Cellular Networks

Intelligent reflecting surfaces (IRSs) are a promising technology for future-generation wireless networks by extending coverage region to blind spots and increasing mmWave propagation paths in non-line of sight environments. User association (UA) in dense millimeter wave (mmWave) networks is vital to characterizing connections among base stations (BSs) and mobile users for load balancing, interference management, and maximizing network utility. However, it has yet to be examined thoroughly in a multi-IRS-aided network. This paper presents a new UA scheme that takes cell interference into account for a multi-cell mmWave cellular network aided with multiple IRSs. We formulate a network spectral efficiency maximization problem by jointly optimizing active beamforming (AB) at the BSs, passive beamforming (PB) at the IRSs, and user-BS association with consideration of the impact of IRSs. We then propose a computationally efficient iterative algorithm based on alternating optimization (AO) to resolve this intractable mixed-integer non-convex problem. A fractional programming technique is used to optimize active beamforming at the BSs and passive beamforming at the IRSs, and a penalization method combined with successive convex programming is applied for UA optimization, which is shown to reach the optimal solution. Simulation results show significant performance improvements obtained by the proposed algorithm, providing higher spectral efficiency compared to several benchmark algorithms, while having a low computational complexity.

acute. In part due to the ability of millimeter waves (mmWave) to provide gigabit data rates and a vast range of bandwidths available between 30 and 300 GHz, millimeter wave communications are considered as part of next-generation cellular networks (5G and beyond) [1]. However, there is a fundamental difference between mmWave and microwave communications where most current communications systems operate. The path loss attenuation, lack of line of sight (LoS) propagation, and sensitivity to blockage all intensify and require careful design in mmWave communications.
MmWave network rollout can be costly because of the dense network requirement to overcome high path loss and blockage sensitivity, together with the MIMO hardware for beamforming. By reconfiguring the propagation environment of wireless communication, the use of intelligent reflecting surfaces (IRS) in wireless communications represents a revolutionary approach to achieving better spectrum and energy efficiency, especially in mmWave networks [2], [3]. IRSs are composed of many low-cost passive arrays that adjust the phase shifters independently to reflect the incident signal. As a result of this controllable reflected signals, the desired wireless transmission can be achieved. With proper phase shift adjustments (referred to as passive beamforming), both the reflected and desired signals can be added constructively and coherently at the receiver, resulting in an improvement in SINR. By taking advantage of the IRS' superior ability to shape wireless propagation environments, wireless networks can be designed with new degrees of freedom (DoFs), resulting in significant improvements in performance. As a passive array, IRS operation differs significantly from that of other active arrays, such as relays, in that IRS requires no additional power and can operate in the full-duplex mode. An IRS located in the LoS with a BS (which is usually the case) improves the signal strength in its vicinity and extends the coverage range. Research efforts have been devoted to designing IRS-assisted wireless communication systems based on these appealing advantages, utilizing a variety of metrics, including maximizing energy efficiency [4], [5], spectral efficiency [6], and the network sum rate [7], [8], [9]. IRS can also be applied in combination with other emerging technologies, such as integrated sensing and communications (ISAC) [10].
User association (UA) has a significant influence on the network performance and interference mitigation.  a mobile user (MU) is assigned to the base station (BS) which provides the maximum signal-to-interference-plus-noise ratio (SINR). As a result, a BS might serve many MUs and suffer overloading, while the other BSs might serve only a few MUs and might be underloaded. UA has been studied in different scenarios, such as massive multiple-input multiple-output (MIMO) systems [11], heterogeneous networks [12], and mmWave Wi-Fi networks [13], and recently in mmWave cellular networks with load balancing [14]. However, UA in a multi-cell mmWave network aided with multiple IRSs has not been fully examined in the literature, considering the effects of IRSs' passive beamforming on the effective channel between each MU and its associated BS.

A. Background and Related Work
Leveraging IRS in a multi-cell network severely affects channel reconfiguration, such as the power allocation and UA algorithms. The effects of IRSs require a new design for an effective active and passive beamforming and UA algorithm to improve the performance of the mmWave cellular network. This joint optimization of active and passive beamforming and UA is challenging and complex. In the context of IRS-aided multicellular systems, UA is still a relatively new topic that has been studied only in a few papers [15], [16], [17], [18], [19], [20], [21], [22]. In Table I, we summarize the representative works on IRS in multi-cell networks based on their considered system setup, design objective, optimization parameters, and frequency band. A sum rate maximization under quality-of-service (QoS) constraints guaranteeing the minimum rate was conducted in [15] using zero-forcing (ZF) receivers. Since ZF receivers suppress the interference, only the signal-to-noise ratio (SNR) is used for obtaining the users' achievable rate. The problem of optimizing the reflective phase shifters at the IRS, the BS's transmit power, decoding order, and subchannel assignment in non-orthogonal multiple access (NOMA) networks was studied in [16]. However, the algorithm requires accurate channel state information (CSI) and has high computational complexity using the many-to-one matching method. The SINR in a multi-IRS and multi-BS network based on the Rayleigh fading model was derived and used in an alternating optimization (AO) to optimize the IRS-user association and BS power allocation jointly [17].
Typically, the IRS and the serving BS have a LoS channel due to the IRS deployment location. Therefore, channel distributions may not always follow Rayleigh fading models. It is also critical to investigate the BS-user association optimization based on arbitrary channel response models in order to evaluate the upperbound performance. Furthermore, optimization-based methods have not been used to examine the impact of the IRS on multi-cell networks and UA. The IRSs' reconfiguration of channels can affect the UA algorithm in a multi-cell environment. Hence, multi-cell networks require joint IRS phase optimization and UA to enhance the performance [18], [19], [20]. Signal attenuation due to the BS-MU distance is one of the most crucial issues in mmWave communications. In our initial work [21], UA in an IRS-assisted mmWave cellular network considering the impact of IRS on spectral efficiency improvement is investigated. By optimizing the UA using a matching game, we balanced the BS loads and maximized network spectral efficiency. We have also considered phase optimization and joint UA for IRS-assisted multi-cell networks in [22].

B. Our Contributions
Most existing works in IRS-assisted communications only consider single BS scenarios, which is not practical to implement in cellular communications with multiple BSs and cells. The system design for IRS-assisted mmWave multi-cell network is a problem worthy of further investigation. Furthermore, the user association problem becomes more challenging and complicated by the IRS deployment at each cell since the cascaded channels (BS-IRS-MU) must be considered in addition to the direct channels. This paper proposes an interference-aware UA scheme for mmWave multiple-input single-output (MISO) cellular networks, in which IRSs are used for coverage enhancement, particularly at the cell edge, and reduce the vulnerability of mmWave to N-LoS paths. We study the joint design of passive beamforming at the IRSs, active beamforming at the BSs, and the UA optimization problem in IRS-assisted mmWave MISO cellular networks. Although some preliminary designs and results have been published in a conference version [22], this paper takes into account both active and passive beamforming in mmWave cellular networks as a means of improving system performance and provides a thorough and detailed analysis of the effectiveness of IRS to improving the system spectral efficiency and cell edge coverage.
This paper makes the following main contributions: r Our study presents an IRS-assisted multi-cell mmWave cellular network, in which an IRS is assigned to each BS in order to facilitate downlink MISO transmission. With respect to the proposed system, we jointly consider active and passive beamforming and UA in order to formulate the problem of sum rate maximization.
r Due to the nonconvexity and complexity of the problem, we present an iterative algorithm for UA, active beamforming, and passive beamforming (IUA/PB) based on AO techniques. The original problem is decomposed into three non-convex subproblems, including UA optimization, active beamforming optimization at the BSs, and reflection phase optimization (passive beamforming) at the IRSs. Each individual non-convex subproblem, is challenging, and we propose two effective approaches to solve them.
r The UA optimization problem is a mixed integer nonlinear programming (MINLP), and due to its non-convexity structure with the presence of the integer variables, is known to be NP-hard. First, the binary UA variables are relaxed using a penalty-based method, and the objective function is approximated by leveraging a first-order lower bound. Next, we design an iterative algorithm based on successive convex programming, which converges quickly to the optimal point due to its low complexity structure.
r To optimize active and passive beamforming, we first deduce the objective function based on the Lagrangian dual transform (LDT). Then by applying the fractional programming method and reformulating the problem as a quadratically constrained quadratic program (QCQP), we design an algorithm to solve the problems effectively.
r Finally, the performance of the proposed joint algorithm is analyzed and compared to different benchmarking algorithms. We show through simulations that the proposed algorithm provides significant improvements in terms of spectral efficiency (SE). The proposed algorithm, compared to benchmarks and a system without IRS, provides a higher level of coverage and lower probability of outages, especially for the cell edge users.

C. Paper Organization and Notation
The rest of the paper is organized as follows. The next section describes the system model and formulates the optimization problem. The proposed IUA/PB algorithm on the basis of AO is designed to solve the problem in Section III. The simulation results of the proposed algorithm are shown in Section IV, and the conclusion is drawn in Section V. Throughout this paper, vectors and matrices are represented by bold lowercase letters and bold uppercase letters, respectively. Superscript (·) * denotes the conjugate transpose. The distribution of a complex Gaussian vector x with mean μ and covariance matrix Q is denoted by x ∼ CN (μ, Q). The expected value operator is denoted by E[·] throughout the text.

A. Mmwave Channel Model
There is a considerable difference in mmWave channels' characteristics and typical microwave channels, which is independent and identically distributed (i.i.d) channel with rich scattering. As is shown in this paper, the mmWave channel model is on the basis of the cluster channel model proposed in [23], which consists of C clusters, where each cluster involves L rays as follows where γ c is the gain of c th cluster and azimuth and elevation angles of arrival and departure are indicated by φ R c,l , ψ R c,l , φ T c,l , ψ T c,l , respectively. a(φ, ψ) represents the uniform planar array (UPA)'s response vector enables 3D beamforming in both azimuth and elevation directions.
Furthermore, two link states, LoS and N-LoS, are considered with probability obtained from the measurements in [24] as where d is the transmitter-receiver distance in meter, d BP is the breakpoint distance shows that there is no longer a %100 probability of LoS, and η is a decay parameter. Based on the measurements, d BP = 27 m and η = 71 m are considered. Moreover, the path loss model for both LoS and N-LoS links can be expressed as where d 0 is the reference distance, λ is the wavelength, n is the path loss exponent, and X σ SF indicates the shadowing lognormal random variable with standard deviation as σ SF . Based on the LoS or N-LoS path, these parameters vary at 28 GHz can be expressed as: n LoS = 1.7, n N −LoS = 3, σ SF,LoS = 3.6 dB, and σ SF,N −LoS = 9.7 dB.

B. IRS-Assisted Mmwave Cellular Network
As illustrated in Fig. 1, we present an IRS-assisted mmWave cellular network, in which several single-antenna MUs are served by several multi-antenna BSs with the assistance of IRSs. Indoor applications with dense users, such as exhibition centers, stadiums, and shopping malls can benefit from such a cellular network [25]. The coverage area of each cellular network is relatively small due to the vulnerability of the BS-MU path to blockage in mmWave communications. Therefore, assigning an IRS to each cell extends the coverage region without consuming any additional power. In this scenario, an IRS is mounted in the LoS path of each BS to reflect the BS's signal and provide an additional path to compensate N-LoS vulnerability in mmWave communications. Hence, the signal is received from two links: direct path (BS-MU) and the reflected path (IRS-MU). Compared to conventional cellular networks without IRS, SINR can be significantly improved when an IRS is applied to each BS.
A downlink mmWave MISO cellular network with J M − antenna BSs and K MUs is assumed. Additionally, an IRS with N phase shifters is assigned to each BS. The sets of BSs, MUs, and IRS phase shifters are denoted by J = {1, . . ., J}, K = {1, . . ., K}, and N = {1, . . ., N}, respectively. Also define Q j as the Activation Set, which provides the set of active MUs served by BS j such that Q j ⊂ K and |Q j | ≤ K. The transmitted signal from BS j is as where s k is a zero mean and unit variance i.i.d random variable denoted as data for MU k and f k,j ∈ C M is a beamformer vector optimized for each MU k by BS j. The power constraint at BS j is expressed as where P is the transmit power at each BS. A multiuser mmWave cellular network with J BSs is shown in Fig. 1, where BS j, j ∈ J, is associated to IRS j with N reflecting passive arrays. The IRSs are capable of dynamically adjusting the phase shifts of their reflecting passive arrays and installed to facilitate communication between the BS and the MU. At MU k, the received signal is expressed as follows where h k,j ∈ C M shows the channel vector between BS j and MU k and z k is the white Gaussian noise at MU k such that . . , θ j,N ] T as a phase-shift vector with reflection coefficients (θ j,n ), and Θ j = diag(θ j,1 , θ j,2 , . . . , θ j,N ) as the diagonal phase-shift matrix (n ∈ N, j ∈ J). Therefore, h k,j is as In terms of interference, the cellular network structure and the UA have a significant impact. Therefore, besides considering fast channel variations in the mmWave [26], interference should be taken into account in UA. The received signal at MU k can be expressed as where the first term represents the desired received signal from BS j. Interference caused by transmitted signals from BS j designed for the other MUs associated with BS j is presented as the second term, interference signals sent by other BSs (z = j) to serve MUs connected with them is represented as the third term, and the last term represents noise received at MU k. The activation sets Q j and Q z significantly influence the interference terms, as shown in (9). As a result, network interference is determined by the network UA. Assuming that MU k is connected to BS j during a transmission, regarding the received SINR at MU k, the instantaneous rate is determined as follows where a k,j is defined as the binary association variable In order for the association variables to be meaningful, they must satisfy the following conditions C1: C2: As reflected in constraint C1, each MU must be served by one and only one BS at a given time. In addition, constraint C2 represents resource allocation across the network in such a manner as to avoid queuing and congestion at the BSs by preventing the number of MUs assigned to each BS from exceeding its number of antennas. As illustrated in (10), the association variables affect instantaneous rate. Thus, the instantaneous rate for MU k is expressed as

C. Problem Formulation
By defining the network's instantaneous rate vector as r (r 1 , r 2 , . . . , r k ), we find the optimized association variables, beamforming matrices, and phase shift matrices that maximize the network's sum rate. Denote a = {a k,j |∀j ∈ J, ∀k ∈ K}, f = {f k,j |∀j ∈ J, ∀k ∈ K}, Θ = {Θ j |∀j ∈ J}, and θ = {θ j |∀j ∈ J} as the set of the association variables for all MUs to the BSs, set of all active beamforming matrices of the BSs, set of all reflective phase matrices (passive beamforming matrices) of the IRSs, and set of all reflective phase vectors of the IRSs, respectively. The optimization problem framework is presented as As described earlier in (12) and (13), C1 and C2 ensure that all the MUs are served by their associated BSs with no need for scheduling. Also, C3 is due to the binary nature of the UA variables, C4 relates to the transmit power constraint, and C5 represents the IRS phase shift bounds. There is a highly nonconvex MINLP problem in (15), known to be NP-hard. MINLP problems are challenging to solve because of their combinatorial structure and the possibility of multiple local maxima in the search space. By using an AO method, the following section will enhance the solution iteratively to the optimization problem in (15).

III. PROPOSED IUA/PB ALGORITHM
Due to the non-convexity of the objective function and the constraints, it is not easy to find the optimum global point of the optimization problem in (15). Therefore, based on the AO technique, a tractable algorithm is proposed to iteratively and separately solve a, f , Θ. In particular, for given f and Θ, a is optimized by applying successive convex approximation. A fractional programming technique is then used to optimize f with a fixed a and Θ using the Lagrangian dual transform to decouple the objective function. Next, we solve Θ by utilizing the fractional programming method. Iterating this process of finding a solution to the problem improves the sum rate at each iteration since, on the feasible set, (15) is upper-bounded. The objective eventually converges to the approximate optimum value. Finally, the IUA/PB algorithm is proposed, and its complexity is analyzed. Note that we use the instantaneous channel knowledge in the following analysis and simulation to develop an upper performance bound for the MU-IRS-BS association. However, the proposed algorithms are not specific to instantaneous channel knowledge and can be applied using other types of CSI, including long-term statistical CSI or a CIS mixture as discussed in more detail in the complexity analysis section.

A. User Association Optimization
This subsection focuses on problem (15) to optimize the association variables a. Considering fixed f and Θ, our focus is on solving the assignment problem Since problem (16) is an integer problem and hence is combinatoric in nature, solving it optimally will require an exponential complexity. Here we design a customized, low complexity algorithm involving two steps to solve this problem. The first step is relaxing the integer constraint by using regularization, which results in a continuous but non-convex optimization problem. The second step is to apply a series of first-order convex lower bound (Lemma 1) to the non-convex objective function to create a sequence of linear optimization problems which can be solved very efficiently with low complexity. We show (in Theorem 1) that the solution of this sequence of linear problems converges to the optimal solution of the regularized problem. By setting a suitable value for the regularization parameter, the regularized problem will yield a solution near the optimal integer association variable. These two steps are described next.

1) Integer Relaxation via Regularization:
There is a challenge in solving the optimization problem in (16) analytically because it is not convex nor concave. To achieve a solution, we propose a heuristic algorithm. First, the binary nature of UA variables (a) is investigated. It can be proved that a k,j ∈ {0, 1}, is equivalent to a 2 k,j = a k,j . Moreover, it holds that in range a k,j ∈ [0, 1], a 2 k,j ≤ a k,j . A penalty term is added to the utility function in (16) enforcing a 2 k,j = a k,j . Hence, the binary integer constraints will be relaxed to a k,j ∈ [0, 1]. Therefore, we have the following optimization problem to solve which is equivalent to (16) when λ >> 0 [27]. Note that the constant λ signifies the importance of the recovering binary variables for a k,j over maximizing the utility function. Also, since the term (a 2 k,j − a k,j ) is always non-positive, the regulation term defined above can be considered as the degree of satisfaction of the binary constraints a k,j ∈ {0, 1} when (17) is solved for λ = ∞ in practice. In other words, in the objective function of (17), the regulation term is added to achieve the optimization results to a binary solution for a k,j . Hence, for λ >> 0, the relaxed problem and the main problem obtain equivalent results [27], [28], [29].
2) Successive Convex Approximation Method for the Regularized Problem: For the purpose of solving the optimization problem in (17), the successive convex approximation (SCA) technique is considered. Due to the affine nature of the constraints C1, C2, and C6, only a lower bound for the objective function is required for applying the SCA method. Note that for a convex function f ( The lower bound of f (x) is affine and can be used to approximate f (x) in the maximization problem using the SCA method. Note that the utility function in (16), U (r), can be expressed as where n k,j , d k,j , and b k,j are as follows Since f (x, y) = x ln(1 + x y ) are jointly convex in x and y, U (r) will be convex (the proof of the convexity is derived in Appendix A). Therefore, the lower bound of U (r) can be first-order approximated at point a (î) . Furthermore, the regulation term in (17) is the sum of quadratic functions in terms of the association variables, which is also convex.
Lemma 1: The first order lower bound of F 1 at a given point a (î) is shown in (20) Proof: The proof is given in Appendix B. A global lower bound maximization for (17) can be expressed as the following convex problem Since problem P1.1 is a linear optimization, it can be solved very efficiently with low complexity. Using linear programming methods, such as the interior point, it is easy to solve the optimization problem in (21). Considering (17), (21),F is the global lower bound of F 1 (a), i.e.
It is, therefore, possible to replace the nonconvex problem in (17) with a sequence of global lower bound linear maximization problem in (21) as follows: first, initialize it from a feasible point a (0) of the problem (17), after that, the optimization problem in (21) is solved iteratively to generate a sequence {a (î) },î = 1, 2, . . . of feasible and improved points toward the optimal solution of (17). Note that at iterationî, a (î−1) is used as a feasible initial point for solving (21) and obtaining a (î) . Theorem 1: After initializing from a feasible point a (0) , a sequence {a (î) } is obtained by solving the linear optimization problem in (21) iteratively. This sequence provides a series of improved points to the regularized problem in (17) and converges to an optimal KKT point for (17).
Proof: The proof is given in Appendix C.
In order to achieve convergence for the AO framework, a series of convex problems (21) must be solved and repeated. For a given error tolerance ξ > 0, with the initial feasible point a (0) , a finite iterations of the UA algorithm leads to the solution of problem (16).
By applying the proposed UA algorithm, (a 2 k,j − a k,j ) is enforced to be 0 by setting a k,j = 0 or 1, and therefore, the integer constraints are satisfied. Also, the UA algorithm terminates after a limited iterations for a given ξ > 0.

B. Active/passive Beamforming Optimization
This subsection focuses on problem (15) to optimize the active and passive beamforming variables (f , Θ). We solve the beamforming problem by considering fixed a. Here the problem is non-convex, and we apply the Lagrangian Dual Transform technique and Fractional Programming of [30] to arrive at a set of closed-form solutions for the active beamforming f and passive beamforming Θ. These closed-form solutions can then be integrated with an overall iterative algorithm to arrive at a solution to the original problem (P1).

1) Lagrangian Dual Transform:
In order to deal with the logarithm in the objective function of P1, the Lagrangian Dual Transform (LDT) is applied [30]. Therefore, the problem in (15) can be expressed as follows for given association variables a max f ,Θ,π U 1 (f , Θ, π) s.t. C4 : The equivalent utility function is as where π = [π 1 , π 2 , . . . , π K ], and π k is auxiliary variable for decoding SINR k . For given f , Θ, the optimal value of π k can be found as π opt k = SINR k . Thus, by substituting π opt k in (24), the optimization problem is reduced to whereπ k = (1 + π k ). (26) is the multiple-ratio FP summation, and fractional programming techniques can solve the nonconvexity of the problem due to the ratio operation [30]. The following two subsections provide further details on how to solve f by fixing Θ, and to solve Θ by fixing f , respectively.
2) Active BS Beamforming: This subsection investigates how to achieve optimal active beamforming matrix f given fixed Θ for (26). Thus, the optimization problem in (26) becomes as Using quadratic transform, F 2 (f ) is reformulated as [30] F where β = {β k,j |∀k, j} are auxilary variables. The optimization problem in (27) can be reformulated to the following problem over f and β (29) is a biconvex optimization problem. To solve it, one common method is to fix one of f and β, then to solve the convex optimization problem corresponding to the other [31]. The optimal β for given f is obtained by setting ∂F 2a (f ,β) Then, fixing β, the optimal f is where λ j is the dual variable introduced for the power constraints that can be determined by applying the sub-gradient method.

3) Passive IRS Beamforming:
Finally, Θ is optimized in (26) given fixed f . Using h k,j defined in (8), the utility function of (26) is represented as By Since F 3a (θ) is fractional programming, on the basis of the quadratic transform, it can be expressed as follows [30] F 3b (θ, ) where = { k,j |∀k, j} are auxiliary variables. The optimization problem is therefore reformulated as such that and θ are optimized alternatively. By setting ∂F 3b (θ, ) ∂ = 0, it is possible to obtain the optimal for a given θ, i.e., Given the optimal opt , the optimization problem for θ is expressed as where (·) represents the real part of a complex number and Since l k,i,j l * k,i,j for all k, i, j are positive definite matrices, B j is a positive definite matrices and F 3c (θ) is a quadratic concave function of θ. As a result, the problem can only be characterized as QCQP, and its non-convexity can only be attributed to the constraints. As an alternative to non-convex constraints, the following convex quadratic constraints may be substituted as where e n ∈ R N represents an elementary vector involving a one at the n th position. As a result, the convex QCQP can be described as follows The above problem is convex, and it can be reformulated to the dual problem via Lagrange dual decomposition (LDD) as where λ = {λ j,n |∀j, n}, λ j,n is the dual variable for the constraint θ j * e n e n * θ j ≤ 1 and L(θ, λ) denotes the dual objective function given by λ) is a concave function of θ satisfying the Slater's condition; thus, the duality gap is indeed zero [32]. Therefore, the optimal θ for a given λ can be obtained by setting ∂L ∂θ = 0 as Proposing a closed-form solution for dual variables (λ) is infeasible. However, it can be determined by applying the sub-gradient method to update λ according to the constraints in (40).

C. Computational Complexity Analysis
Algorithm 1 summarizes our proposed IUA/PB algorithm for the initial problem (15). In particular, the association variable and formulation can work with any time scale, as long as the SINR in the sum rate used for the association problem (P1.1) is updated at that time scale. For the proposed algorithms, different optimization parts can work at different CSI time scales, for example, UA (problem P1.1) with statistical CSI while beamforming and IRS phase shifts (problem P1.2, P1.3) with instantaneous CSI. As stated before, our approach to solving a, f , Θ consists of continuous and alternative iterations until a stable optimal point is reached. Algorithm 1 is guaranteed to converge since each optimization increases the sum rate value after each iteration, and the objective is an upper bound over the feasible set of (15).
In general, the complexity of the proposed algorithm varies with the number of iterations in the outermost layer alternation and with the complexity required to solve each subproblem at each iteration. For the UA subproblem, the computational complexity of solving (21) is polynomial in the number of variables and constraints. The (21) is an optimization problem with a (KJ) real-valued scalar decision variables, a linear objective, and b (K + J + KJ) linear constraints. Therefore, the complexity required to solve (21) is O ((1 + a + b)

IV. SIMULATION RESULTS
An evaluation of the effects of IRS and the performance analysis of the proposed IUA/PB algorithm are presented in this section in a simulated downlink mmWave cellular network that employs J BSs and K MUs operating at 28 GHz. Based on the explanation in Section II-A, mmWave channels are generated such that there are 5 clusters in each channel, with 10 rays per cluster. There are 8 × 8 UPA antennas installed at each BS, 5 × 5 arrays installed at each IRS (the number of passive elements in the IRS needs to be adequate to achieve an effective indirect path), and a single antenna at each MU. BSs transmit with equal power P , and the power spectral density of noise is −174 dBm/Hz. While the BSs' location are fixed and known, the IRSs are located in random locations surrounding the BSs, and the MUs are randomly distributed within a region of 100m × 100 m. It is necessary to consider a LoS path between each BS and its corresponding IRS in order for the IRS to be effectively utilized. Furthermore, D j is the maximum number of allowed active mobile users at each BS. If the number of MUs associated with a BS is more than D j , the BS suffers overloading, and congestion occurs.
To evaluate our proposed IUA/PB algorithm, we compare it with four algorithms under various scenarios: r DA-MG+PPB+MRT: Deferred Acceptance-Matching Game (DA-MG) method [34], [35] for user association, the proposed passive beamforming (PPB) at the IRSs, and maximum ratio transmission (MRT) at the BSs are applied.
r Max-SINR+PPB+PAB: The conventional max-SINR scheme [14] is utilized for UA. Also, PPB and proposed active beamforming (PAB) are applied at the IRSs and BSs, respectively. r PUA+RPB+PAB: The proposed user association (PUA) and PAB are considered, while the phases of the IRS arrays are assigned randomly to achieve random passive beamforming (RPB).
r PUA+No-IRS+PAB: While PUA and PAB are applied, no IRS is considered in BS/MU communications. In this simulation, we consider 5 mmWave BSs and 25 MUs, while an IRS is associated with each BS for improving the communications link. Also, it is assumed that each BS can serve at most 5 MUs simultaneously.

A. Spectral and Energy Efficiency
For the different UA schemes with and without IRS assistance in the network, Fig. 2 shows the network spectral efficiency (SE) based on the sum-rate utility function U (r) given in (15). The figure illustrates that our proposed IUA/PB algorithm outperforms and the network interference is adapted accordingly. Furthermore, since joint UA, active beamforming, and passive beamforming optimization are performed simultaneously in the proposed scheme, it achieves higher SE in comparison with the other schemes. Specifically, when P = 10 dBm, based on our proposed algorithm, we can achieve 109% sum rate in comparison to DA-MG+PPB+MRT and 116% sum rate in comparison to Max-SINR+PPB+PAB. Also, the impact of the IRS in the mmWave cellular network is illustrated in Fig. 2. IRS is mostly useful at mmWave frequencies, where there may be a sparse channel and no LoS in the BS-MU link. Thus, an additional propagation path through IRS is required, despite the fact that this path is weak as a result of channel attenuation at this frequency. The IRS is therefore able to provide a higher sum rate when applied to cellular networks.  Similarly, Fig. 3 illustrates the energy efficiency (EE) performance trend under the same settings as Fig. 2. Note that the EE of the system is expressed as the network spectral efficiency over the network power consumption as where P cir is the total circuit power in the cellular network considered as 20 dBm. As shown in the figure, our proposed IUA/PB algorithm can provide the best EE, where the EE of the IUA/PB algorithm is 9% greater than that of DA-MG+PPB+MRT and 16% higher than that of Max-SINR+PPB+PAB, respectively. This is due to the joint optimization gain of UA and active/passive beamforming. Moreover, all of these algorithms have an increased EE with increasing transmit power in the low transmit power range, similar to the trend of Fig. 2. EE performance is primarily influenced by the SE rather than the transmit power in the low region. However, after increasing from a specific transmit power value, we can  Fig. 4(a)) and at each IRS ( Fig. 4(b)).
foresee that since the excess BS transmit power is not utilized in the existing mmWave system, the EE performance is negatively affected by increased transmit power. As a result, our algorithm offers significant EE gains within a range of acceptable sum rates.

B. Impact of BS and IRS Antenna Sizes
The impact of array size on the network spectral efficiency is shown in Fig. 4, which illustrates the network spectral efficiency versus array size at each BS, and each IRS. In Fig. 4(a), the number of antennas at each MU and IRS are fixed (N = 5 × 5), and the number of BSs' antenna is increased. It is shown that by increasing the antenna array size from 3 × 3 UPA to 10 × 10 UPA, the network spectral efficiency of the proposed algorithm improves by 95%. Furthermore, comparing the proposed joint optimization IUA/PB algorithm with the other algorithms, it achieves the highest spectral efficiency. It is because, an increase in M will result in an increase in the number of antennas that can be used for beamforming, thereby improving the efficiency of active beamforming and the spectral efficiency of a network. In addition, the proposed algorithm achieves significant improvements in spectral efficiency due to better UA, active beamforming, and passive beamforming gains, demonstrating the effectiveness of our joint optimization process. In Fig. 4(b), the IRS array size affects the network spectral efficiency while the array size at BS is fixed (M = 8 × 8). The greater the number of IRS arrays N , the greater the network's spectral efficiency. In general, as the number of reflective arrays increases, more signals can be reflected, and therefore more effective passive beamforming can be achieved. The IRS may also have a significant effect when there is a sufficient number of phase shifters. All the results confirm that increasing the array sizes at either BS or IRS significantly impacts the network spectral efficiency. Also note that since random phase beamforming is applied in the PUA+RPB+PAB method, the reflected links will not be added constructively to the desired signals. Furthermore, the wireless propagation environment created by random IRS phases might be shaped differently than the desired one. These non-constructive reflected links are more apparent when the number of passive elements is not large enough (such as 4 × 4) to provide a significant reflected link. Therefore, applying IRS without optimizing its phases might be useless, and the system performance might be the same without IRS. Fig. 5 illustrates how the number of MUs (K) impacts the SE of the network. In general, there is an apparent ascending trend among all five schemes when K increases within a specific range due to providing higher network resource utilization. Despite this, power resources in BSs limit the SE's improvement indefinitely. In addition, our proposed algorithm has a higher performance than other algorithms by taking advantage of the optimal gain brought about by active beamforming at the BS, passive beamforming at IRS, and UA algorithms. As shown in Fig. 5, when K = 50, based on our proposed algorithm, we can achieve a 6%, 13%, 21%, and 45% improvement in performance compared with DA-MG+PPB+MRT, Max-SINR+PPB+PAB, PUA+RPB+PAB, and PUA+No-IRS+PAB, respectively.

C. Achievable Rate Distributions
Cell edge users typically experience the lowest transmission rates. The probability density function (PDF) and the cumulative distribution function (CDF) of the users' data rate in the mmWave cellular network for different schemes are depicted in   6. It can be inferred that the probability of the users having a low rate becomes small by applying the proposed IUA/PB algorithm in comparison with the other algorithms. Due to the fact that most users located in the cell-edge would belong to the low-rate region, these results confirm that using the proposed IUA/PB algorithm increases the data rate for cell-edge users. As shown in these figures, the low-rate region occurs with very low probability (lowest among all schemes) because of the improved link quality provided by the IRS and the proposed IUA/PB algorithm. Also, using the IRSs and our proposed algorithm increases the high transmission rates significantly, evident by the probability density function reaching a much higher rate range. This result shows that our proposed algorithm helps to both improve the cell edge users, who typically suffer low rates and increase the probability of a user achieving a high transmission rate.

D. Algorithm Convergence
Finally, Fig. 7 shows the convergence of Algorithm 1 for the proposed IUA/PB algorithm with 5 BSs, K MUs, and N passive elements (K = 25, 50, 75; N = 5 × 5, 8 × 8). As shown in the figure, the proposed algorithm is fast to converge and requires  only a maximum of five iterations to achieve convergence. There are only three convex problems with polynomial complexity to be solved at each iteration of the proposed IUA/PB algorithm. An average is calculated over 500 channel realizations, and the elapsed time of the proposed algorithm is calculated. As expected, based on our complexity analysis in Section III-C, the actual elapsed time is raised by increasing the number of MUs and IRS passive arrays due to the complexity increments. Fig. 8(a) shows the convergence of the objective function (21) in the inner loop (î) for the user association optimization sub-problem in (16). By using the regularization method to transform (16) into (17), the regularization parameter needs to be large enough to enforce the integer constraint. However, at the same time, a larger λ will result in a worse original objective function, and Fig. 8(a) illustrates this trade-off. Note that due to the penalty term in (17), which is non-positive, the objective function at the first iteration gets negative. By iteratively enforcing the regulation term to zero, the algorithm converges to a positive sum-rate value. Fig. 8(b) shows after some iterations, the optimal association variables become very close to either 0 or 1 to force the regulation term toward zero and make the integer constraints valid. Also, there is a trade-off between having a large λ to enforce the association variable to be an integer with fewer iterations for convergence and having a small λ with many iterations for convergence and achieving an exact optimal point. To make a compromise between them, in our simulations, λ is chosen to be 5. Fig. 8(b) displays a scatter plot of the distribution of one of the association variables, which ultimately converges to 1 based on different λ. As the figure shows, the convergence rate depends on the value of λ. Larger values of λ result in faster convergence with fewer iterations.

V. CONCLUSION
To realize cost-effective and high-coverage mmWave communication, we considered an IRS-assisted mmWave cellular network. Considering the IRS effect on UA, a spectral efficiency maximization problem is formulated by jointly optimizing UA with load balancing constraints, active beamforming at the BSs under transmit power constraints, and passive beamforming at the IRSs under phase constraints. Since the formulated problem is non-convex and mathematically intractable, by combining successive convex programming with the penalty method and fractional programming, an iterative method on the basis of alternative optimization has been developed. The proposed joint algorithm improved the objective spectral efficiency in each iteration and was shown to converge. The feasibility and effectiveness of our proposed algorithm were demonstrated through extensive simulations. Specifically, as compared with the benchmark, the proposed algorithm is capable of providing significantly higher SE, while having low computation complexity by achieving convergence within only a few iterations. These results also showed that the use of IRSs helped improve coverage probability for low-rate users by significantly boosting their spectral efficiency. APPENDIX A PROOF OF THE CONVEXITY OF f (x, y) = x ln(1 + x y ) Calculating the 2 × 2 Hessian matrix elements of the function x ln(1 + x y ) yields The determinant of Hessian matrix is det(Hessian) = 0.
Therefore, one of the Hessian matrix eigenvalues is zero, and the other one is always positive since the summation of all diagonal elements of the Hessian matrix is always positive for x, y ≥ 0.

APPENDIX B THE LOWER BOUND APPROXIMATION OF F 1
Since U (r) in (18) is convex, its lower bound can be first-order approximated at point a (î) as follows [36] k∈K j∈J Also, (a 2 k,j − a k,j ) is a convex quadratic function, and its lower bound can be approximated as therefore, considering (50) and (51) at the point a (î) , F 1 (a) can be approximated as given in (20).