{'fraction_non_alphanumeric': 0.02989217994794894, 'fraction_numerical': 0.011426446895526087, 'mean_word_length': 4.406861431251675, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimized Power Control for Massive MIMO with Underlaid D2D Communications

In this paper, we consider device-to-device (D2D) communication that is underlaid in a multi-cell massive multiple-input multiple-output (MIMO) system and propose a new framework for power control and pilot allocation. In this scheme, the cellular users (CUs) in each cell get orthogonal pilots which are reused with reuse factor one across cells, while all the D2D pairs share another set of orthogonal pilots. We derive a closed-form capacity lower bound for the CUs with different receive processing schemes. In addition, we derive a capacity lower bound for the D2D receivers and a closed-form approximation of it. We provide power control algorithms to maximize the minimum spectral efficiency (SE) and maximize the product of the signal-to-interference-plus-noise ratios in the network. Different from prior works, in our proposed power control schemes, we consider joint pilot and data transmission optimization. Finally, we provide a numerical evaluation where we compare our proposed power control schemes with the maximum transmit power case and the case of conventional multi-cell massive MIMO without D2D communication. Based on the provided results, we conclude that our proposed scheme increases the sum SE of multi-cell massive MIMO networks.


I. INTRODUCTION
Device-to-device (D2D) underlay communication and massive multiple-input multiple-output (MIMO) are two new promising technologies in wireless communication that will appear in 5G networks [2]. Although these two technologies show significant improvements in the energy and spectral efficiency of the network, their combination can even enhance the network performance further. D2D underlay communication enhances the spectrum utilization by reusing the cellular resources for direct communication between D2D pairs when the transmitter and receiver are closely located. It provides benefits such as cellular traffic offloading for cellular networks. This gain is due to the direct communication between the transmitter and receiver of D2D pairs instead of sending the data through cellular base stations (BS); higher data rate and lower transmission power between D2D users due to the short-range communication [3]- [5]. Hence, D2D underlay communication increases the spectral and energy efficiency of the cellular networks. These benefits come at the cost of causing extra interference to the cellular users (CUs) [6]- [8].
A preliminary version of this work was presented at the 22nd international ITG workshop on smart antennas (WSA 2018) [1]. This paper was supported by the European Union's Horizon 2020 research and innovation programme under grant agreement No 641985 (5Gwireless).
Massive MIMO is one of the most significant technologies in 5G, as it offers large improvements in spectral and energy efficiency of cellular networks. Because it utilizes a large number of antennas at each BS, it offers multiplexing gains and spatial interference suppression for CUs [9]- [11]. By mitigating the extra interference that the D2D communication causes to the cellular network, we can potentially enhance the mentioned benefits of D2D underlay communication to the cellular network. Massive MIMO seems to address the shortcoming of D2D underlay communication. Hence, the combination of these two technologies has received considerable attention. The idea of combining D2D underlay and massive MIMO technologies has been investigated in both single-cell and multi-cell setups with various objectives, such as data power optimization, proposing new pilot allocation scheme, etc. Although, downlink (DL) D2D underlay massive MIMO has been considered in the literature [12], [13], most of these works focus on the uplink (UL) D2D underlay setup. This is due to the fact that interference suppression is easier to implement in the UL at the massive MIMO BS.
Most prior works in D2D underlaying UL of massive MIMO systems has been investigated mainly for single-cell setups [14]- [18]. In these papers, the transmitter and receiver of a D2D pair are always physically located in the same cell. This simplifies the resource allocation since each D2D pair shares resources with that specific cell. These papers consider problems such as pilot reuse schemes for D2D underlay [14]- [18], data power control for D2D users only [15]- [17], and data power control for both D2D and cellular users [18]. The multi-cell scenario has been studied in a few papers such as [19], [20].
Lin et al. [19] investigate the interplay between massive MIMO and underlaid D2D for UL data transmission in a multi-cell massive MIMO setup where the number of D2D pairs in each cell follows a Poisson distribution. A detailed analysis is provided for the case of perfect channel state information (CSI), in which they study the asymptotic and non-asymptotic spectral efficiency (SE) of CUs and D2D pairs. They also perform channel estimation based on an orthogonal training scheme that allocates orthogonal pilots to CUs and a limited number of D2D pairs nearest to the location of the BS in their cell. In this case, they study the asymptotic and non-asymptotic SE of CUs only. It is assumed that all users transmit with a predefined transmit power. In [20], the authors consider a multi-cell massive MIMO UL D2D underlay model and consider the spatial location of BSs and D2D transmitters follow two independent homogeneous Poisson point process. They derive exact expressions of the SEs for cellular and D2D communication which are not in closed form. Then they apply open-loop power control for interference cancellation in the network. This paper does not consider channel estimation and pilot transmission and they only apply a simple open-loop power control scheme.

A. Contributions of the paper
We consider a multi-cell massive MIMO setup, in contrast to most prior works, with D2D underlaying in the UL. It is assumed that the CUs in a cell have orthogonal pilots that are reused with reuse factor one between cells. In addition, the system contains multiple D2D pairs that are arbitrarily located, do not belong to any cell, and share a network-wide set of orthogonal pilots. Power control is mandatory to achieve good performance in such networks, but this has not been considered in prior works. The existing SE expressions for D2D communication depend on the small-scale fading coefficients and are unsuitable for power control since it should exploit the timefrequency diversity against fading instead of counteracting the instantaneous fading realizations. To overcome these issues, we make the following main contributions: • We derive a new lower bound on the ergodic capacity of the D2D receivers under imperfect CSI. The bound contains an expectation with respect to the small-scale fading. To enable the derivation of tractable power control schemes, we derive a closed-form approximation of the capacity bound and interpret its structure. • We derive closed-form capacity lower bounds for the CUs with either maximum-ratio (MR) and zero-forcing (ZF) processing, which support arbitrary pilot allocation among the D2Ds. • We formulate several data power control problems using our new SE expressions (lower bounds on the capacity) and solve them using convex optimization. We jointly optimize the data power of CUs and D2D pairs to guarantee max-min fairness to the cellular and D2D communications. We also maximize the product of the signalto-interference-plus-noise ratios (SINRs) of the cellular and D2D users, which gives priority to users with good channels. For the first time in this area, we also consider the joint optimization of the data transmission power and pilot transmission power: • We formulate novel joint data and pilot power control problems with MR processing at the BSs and solve both max-min fairness and max product SINR problems using geometric programming. • We propose a successive approximation algorithm to solve the corresponding joint data and pilot power control problem for the case of ZF processing at the BSs. We provide a numerical performance evaluation of our proposed algorithms. The conference version of this work is found in [1], where max-min fairness power control is considered with variable data transmission power only. This paper contains more general and complete results.

II. SYSTEM MODEL
In this paper, we consider a multi-cell massive MIMO system consisting of B BSs, each equipped with an array of M antennas and each BS serves K single-antenna CUs. The system also contains L D2D pairs that are spread over the whole network and do not belong to any specific cell. 1 The network model is illustrated in Fig. 1 with hexagonal cells. Note that this figure is only for illustrative purpose, while our analysis applies to cellular networks with arbitrary shapes of the cells. In this setup, we investigate D2D communication that is underlaid the UL data transmission of the multi-cell massive MIMO system. The wireless channels are varying over time and frequency, which we model by conventional block fading. We define the coherence interval of a channel as the time-frequency block in which the channel is constant and its size is τ c samples (or channel uses) [9,Ch. 2], [10,Ch. 2]. The channels change independently from one block to another according to a stationary ergodic random process. The number of samples per coherence interval is given by τ c = T c B c , according to the Nyquist-Shannon sampling theorem, where T c is the coherence time and B c is the coherence bandwidth. In the proposed setup, indicates the large-scale fading between D2D transmitter l and BS b. In addition, g l,c b,k ∼ CN (0, β l,c b,k ) denotes the channel response between CU k located in cell b and D2D receiver l, while g l,d l ′ ∼ CN(0, β l,d l ′ ) denotes the channel between D2D transmitter l ′ and D2D receiver l. In addition, β l,c b,k and β l,d l ′ refer to the large-scale fading between CU k located in cell b and D2D receiver l and between D2D transmitter l ′ and D2D receiver l, respectively. In our network setup, it is assumed that the BSs and D2D receivers have no prior CSI at the beginning of a coherence interval; hence, channel estimation is carried out in each coherence interval. Therefore, the communication consists of two phases: UL pilot transmission and UL data transmission. To enable channel estimation, the CUs and D2D transmitters transmit pilot sequences of length τ in the pilot transmission phase and the remaining τ c − τ samples are utilized in the data transmission phase. In this setup, we can construct τ mutually orthogonal pilot sequences that are vectors of length τ. It is assumed that we have K orthogonal pilots for CUs that are reused in each cell and we have another set of N orthogonal pilots designated for the L ≥ N D2D pairs. Hence, we have τ = N + K.

A. Uplink data transmission
The received signal during data transmission at BS b is and the received signal at D2D receiver l is where p c b,k and p d l ′ are the transmit powers for data transmission of the CU k in cell b and the D2D transmitter l ′ , respectively. w b ∼ CN (0, I M ) and w l ∼ CN (0, 1) indicate the normalized additive white Gaussian noise at BS b and the lth D2D receiver, respectively. Also, s c b,k and s d l ′ denote zero mean and unit variance data symbols transmitted from CU k in cell b and D2D transmitter l ′ , respectively.

III. ANALYSIS OF SPECTRAL EFFICIENCY
In this section, we derive the SEs achieved when using the communication setup defined in Section II. Pilot contamination arises between CUs due to the pilot reuse between cells and pilot contamination also arises between the D2D pairs that are using the same pilot.

A. Pilot transmission and channel estimation
We denote the matrix of pilot sequences as Φ = [Φ c Φ d ] and it has size τ × (K + N). The matrices Φ c = [φ c 1 , . . . , φ c K ] ∈ C τ×K and Φ d = [φ d 1 , . . . , φ d N ] ∈ C τ×N contain the orthogonal unit-norm pilot sequences assigned for CUs and D2D pairs, respectively. The pilot φ c k is used by CU k in each of the cells.
where p p,c b,k denotes the pilot transmit power used by CU k in cell b and p p,d l is the transmit power of D2D transmitter l for pilot transmission. W b ∈ C M×τ is the normalized additive white Gaussian noise at the BS b which has independent entries having the distribution CN (0, 1). Each BS multiplies the received signal matrix in the pilot transmission phase with each of the pilot signals to despread the signals. Hence, the received pilot signals at BS b from its CUs after despreading by Φ c is [9, Ch. 3], [10, Ch. 3] where W ′′b,c = W b Φ c ∈ C M×K is the normalized additive white Gaussian noise at the BS b with independent entries having the distribution CN (0, 1). The received pilot signal at BS b from all D2D transmitters is the second part of (3) and despreading by multiplication with Φ d yields where W ′b,d = W b,d Φ d ∈ C M×N is the normalized additive white Gaussian noise at the BS b which has independent entries having the distribution CN (0, 1). The minimum meansquared error (MMSE) estimates of the channel vectors from CU k at BS b are [21], [22] The MMSE estimate at BS b of the sum of the channel vec- In addition, the MMSE estimates of the channel vectors from D2D transmitter l ∈ N i at BS b arê The mean square of the channel estimates at the BS are Note that the channel estimateĥ b,d l of D2D transmitter l ∈ N i is a scaled version ofĥ b, N i , which is a property that we will utilize later in the receive processing. The D2D receivers apply the same procedure for channel estimation. Therefore, the received signals at D2D receiver l from the pilot transmission of CU k located in cell b and D2D transmitter l ′ ∈ N i after despreading the signals are where w d l ∼ CN (0, 1) and w c l ∼ CN (0, 1) are the normalized additive white Gaussian noise terms at D2D receiver l. The pilot transmission of all transmitters, i.e., cellular and D2D, are used for channel estimation and the MMSE estimates of the channels from CU k located in cell b and from D2D transmitter l ′ ∈ N i are respectively given bŷ In addition, the mean square of the channel estimates at D2D receiver l are

B. Spectral efficiency with MR processing
In this subsection, we assume that MR processing is applied at all BSs to detect the signals of their own users. The received data signal of user k at the cell b is expressed as (18), as shown at the top of the next page, where the first term is the desired part of the received signal, the second term is cellular intracell interference, the third term is intercell interference, the fourth term is D2D interference, and the last term is noise. We rewrite the received data signal of the CU k at the cell b as where the first term is treated as the desired part of the received signal and the rest of the terms are treated as noise in the signal detection. More precisely, we utilize the use-and-thenforget technique [9, Ch. 3] to lower bound the capacity of each of the CUs, using the capacity bound for a deterministic channel with additive white non-Gaussian noise provided in [9, Sec. 2.3]. We get the lower bound on the capacity of CU k in cell b as where I b,c b,k (MR) is defined as seen in (21) at the top of the next page.
We can write (20) as where I ′b,c b,k (MR) is defined in (23), shown at the top of the next page.

C. Spectral efficiency with zero-forcing processing
In this subsection, we consider ZF processing at all BSs, instead of MR processing, to suppress interference between the users. In addition to suppressing intracell interference, we also suppress interference from the D2D transmitters. For ease of notation, we rewrite (1) in matrix form as Hence, (24) is rewritten as As mentioned previously, the channel estimates at each BS b for D2D transmitters that use the same pilot are parallel vectors. Hence, the collection of N vectorŝ h b,n 1 , . . . ,ĥ b,n N from (7) contains scaled version of all the L channel estimatesĥ b,d 1 , . . . ,ĥ b,d L . Therefore, we construct the M × (K + N) matrix to describe all the channel directions that are relevant for interference mitigation at BS b. We can write this matrix aŝ at the diagonal. By using this matrix, we define the ZF detection matrix at cell b as By applying this decoder matrix in cell b, the received signal becomes as seen in (29) at the top of the next page and the received signal of CU user k at the cell b is provided in (30), that can be seen at the top of the next page, where the first term is the desired signal and the rest is interference plus noise and its variance is provided in (31) can be seen at the top of the next page. Note that this ZF matrix will suppress interference not only between CUs, as in conventional massive MIMO, but also from D2D transmitters. By using the capacity bounding technique in [9, Sec. 2.3], we have the lower bound on the capacity of the CU k in cell b as where I b,c b,k (ZF) is defined in (33) that can be seen at the bottom of the next page.

D. Spectral efficiency of D2D communication
To calculate the SE of the D2D transmissions, we use the MMSE estimates provided in (14) and (15). To detect the desired data from D2D transmitter l, D2D receiver l multiplies the received signal in (2) with the complex conjugate of its channel estimateĝ l,d l . In addition, D2D receiver l uses the channel estimatesĝ l,c b,k ∀b, k andĝ l ′ ,d l ∀l ′ as side information during data detection. We define  (34), is the power of the desired signal. The first term in the denominator is the variance of channel estimation error from the desired D2D transmitter, caused by the gain uncertainty in the detector. The second term is interference from cellular users and the third term includes the interference from other D2D transmitters. This expression is useful for SE computation but it is not in closed-form as it contains an expectation with respect to the small-scale fading coefficients, thus it is not suitable for power control optimization. Hence, to provide a tractable power control algorithm, we compute an approximation of (34) by computing the expectation of the numerator and the denominator of the fraction inside the logarithm. The resulting approximation of the SE of D2D pair l is denoted asR d l and is given by where SINR d l is seen in (36) at the bottom of the next page. Note that in (36) the last expression follows from replacing the mean square of the channel estimates with their actual representation. The tightness of this approximation is investigated in Appendix A.

IV. OPTIMIZATION OF POWER ALLOCATION
In this section, which contains the main contribution of the paper, we investigate different power control schemes for the studied system model. We consider power control where the data transmission powers of CUs and D2D transmitters are selected to maximize the max-min fairness and max product SINR optimization objectives. We also propose a power control scheme for joint pilot and data transmission optimization powers for these two objectives, which is considerable more complicated. In all these cases, we consider both the cases of MR and ZF processing at the BSs, which require different solution algorithms. Note that in all optimization problems provided in this section, P max > 0 denotes the maximum transmit power of the users. In Table I, we summarize the different problems that we solve and the corresponding subsections where the solutions are found. MR processing ZF processing max-min, data power control Section IV-A1 Section IV-A1 max-min, data & pilot power control Section IV-B1 Section IV-C max product SINR, data power control Section IV-A2 Section IV-A2 max product SINR, data & pilot power control Section IV-B2 Section IV-C

A. Data power control
In this subsection, we select the data powers in the system to either maximize the minimum SE in the network (i.e., provide a uniformly high SE among the users) or the product of the SINRs (i.e., the sum of the logarithms of the SINRs, which is a variation of classical sum-SE optimization that guarantees non-zero SE to all users [10, Sec. 7.1]).

1) Max-min fairness:
We will now select the data powers in the system to maximize the minimum SE of the CUs in all cells as well as all the D2D pairs, for both the cases of MR and ZF processing at the BSs. The max-min problem is formulated in epigraph form as where λ indicates the SE level that is guaranteed to all CUs and all D2D pairs and this variable is to be maximized. Note that we use the SE approximation in (35) for the D2D pairs to achieve a tractable problem formulation, while the exact SE will be used later to evaluate the performance. To solve (37),
if the feasibility problem is not satisfied then we fix λ and solve the resulting linear feasibility optimization problem which determines if the given SE level λ is achievable or not. This is a linear program that can be solved efficiently using CVX [23], or some other standard solver for linear programs. To find the optimal λ, we can perform a line search over the interval [0, λ u ], where λ u is an upper bound on the SE level that can be guaranteed to all users. One way to obtain λ u is the compute upper bounds on the SEs achievable by the cellular and D2D users by neglect the interference terms and assuming that everyone transmits at maximum power. The user that achieves the lowest SE under these circumstances limits the SE level that can be guaranteed to all the users. Mathematically, this can be computed as λ u =min l,k,b, in the case of MR processing. With ZF processing, it can be computed as λ u =min l,k,b log 2 1 + P max p d l γ d l,l , We use the bisection line search algorithm for solving problem (37) [24], as described in Algorithm 1. This algorithm can solve problem (37) for both case of MR and ZF processing at the BS, although the optimal solution will be different. Algorithm 1 consists of solving the linear program (37) that has complexity of O((K B + L + 1) 2 (K B + L)), in addition, total number of iteration for convergence of bisection line search is given as ⌈log 2 λ u ǫ ⌉. Therefore, Algorithm 1 has computational complexity order of O(⌈log 2 λ u ǫ ⌉ (K B + L + 1) 2 (K B + L)) [24,Ch. 1]. The exact complexity depends on which solver of linear programs that is used and on the channel realizations that appear in the instance of the problem that is being solved.
2) Max product SINR: The max-min fairness power control approach proposed in the previous section may lead to low overall performance. This is due to the fact that the networkwide performance is limited by the channel condition of the most unfortunate users in the network. When increasing the number of active users in the network, the risk of having an extremely bad channel for some of the active users is high. This is handled in [9] by dropping such users from service before applying the max-min power control. Another potential approach is to perform max product SINR power control instead. To investigate this, we will now change the optimization objective to maximize the product of the SINRs of all users in the network. It offers some level of fairness as well as approximately maximizing the sum SE of the network [10,Sec. 7.1]. The optimization problem is formulated as where for the case of MR processing the SINR b,c b,k in (40b) constraint is given in (22) and for ZF processing we replace it with the SINR expression given in (32). The SINR constraints can be rearranged as geometric constraints, thus the max product SINR with data power control is a geometric program and it can be solved efficiently by using standard convex optimization solvers such as CVX [23]. The interested reader can find the basic terminology and property of geometric programming in Appendix B.

B. Joint pilot and data power control for MR processing
In the previous part, it was assumed that the pilot transmit power of the CUs and D2D transmitters are fixed, e.g., at its maximum value. This is the simplest assumption for pilot transmission and makes perfect sense in a system without pilot contamination. However, in the studied multi-cell scenario with pilot contamination, we can improve the performance by optimizing the pilot transmission powers in addition to data power control. Hence, the goal of this subsection is to perform joint pilot and data power control for the case of MR processing at the BSs with the same optimization objectives as above.
1) Max-min fairness: In this part, similar to Section IV-A1, the optimization objective is to maximize the minimum SE of the CUs in all cells as well as all the D2D pairs. The difference from (37) is that we have extra variables and constraints when we want to optimize the pilot powers as well. The optimization problem is formulated in epigraph form as where λ indicates the SE that is guaranteed to all CUs and all D2D pairs and this variable is to be maximized. By rearranging in the SINR constraints, this optimization problem is in the form of geometric program and as it is mentioned in previous part, this type of optimization problem can be solved efficiently by standard convex optimization solvers such as CVX [23].
2) Max product SINR: For the case of MR processing at the BSs, the optimization problem is formulated as where in the first and second constraints, the SINR expressions are provided from the SE expression given in (22) and (36), respectively. Similar to (40), the SINR constraints can be rearranged as geometric constraints, thus (42) can be solved efficiently as a geometric program.

C. Joint pilot and data power control for ZF processing
In this subsection, we propose an efficient algorithm for pilot and data power control for the case of ZF processing at the BSs. Due to the more complicated structure of the SINRs with ZF, we can only solve the power control problems to local optimality. 2 We focus on max product SINR optimization problem, but it is straightforward to apply the same methodology for the case of max-min fairness optimization problem. The optimization problem of interest is formulated as where the SINR constraint in (43b) is provided in the SE expression given in (32). The optimization constraints (43b) are not posynomial; hence, the optimization problem is not a geometric program (see Appendix B for details). In fact, the SINR constraints are signomial constraints, thus there is no way to reformulate them to make the problem efficiently solvable, with resorting to approximations. The signomial constraints are in the form of fractions of two posynomial functions. To address this problem, we develop an algorithm that finds a local optimum. We propose to first lower bound the posynomials function in the denominator of the signomial constraint with a monomial function. Then, the resulting function is a fraction of a posynomial with the monomial which is also a posynomial function. Consequently, the lower-bounded SINR constraints is a geometric constraint, so that we have obtained an approximation of the original problem that can be solved as a geometric program, using standard solvers such as CVX [23]. By successively updating the approximation, a local optimum is guaranteed to be found. The detailed procedure is described below.
which can be lower bounded by the following monomial functionf (x) where Q j ≥ 0 is a weight. For any fixed and element-wise positive x 0 ,f (x) is the best local monomial approximation of f (x) near x 0 in the sense of first order Taylor approximation if we choose the weight as Proof 1: The proof is provided in [25, Lemma 1]. In the SINR expression of ZF processing in (32), we have (33) in the denominator. In (33), we find which prevent the denominator from being a posynomial function. However, the denominator of f p p,c b ′′ ,k ′′ and g p p,d l , can be lower bounded using Lemma 1 and we have the following upper bound for f p constraint. The second condition ensures the solution from the previous iteration is also feasible for the current iteration while increasing the optimization objective. If the convergence criteria is satisfied the optimal solution of previous iteration is also the solution of current iteration. As the SINR expression is bounded from above and it is a continuous function of p p,c b ′′ ,k ′′ and p p,d l , it ensures that the iterative algorithm converges to a limit point. In addition, the solution of optimization problem with the lower bounded SINR constraint satisfies the Slater's condition. Hence, the limit point of the successive approximation algorithm after convergence is a local optimal KKT point [26], [27]. Each iteration of Algorithm 2 consists of solving optimization problem (43) which has 3(K B + L) variables and the same number of constraints. Hence the order of computational complexity of the algorithm is where N is the required number of iteration for convergence of the algorithm and F is the cost of evaluating the first and second derivatives of the objective and constraint functions in (43) [24,Ch. 1]. Note that the same solution approach can be used for the case of max-min fairness with joint pilot and data power control for ZF processing at the BSs. Therefore, it is omitted to avoid unnecessary repetition in the content.

V. NUMERICAL ANALYSIS
In this section, we provide a numerical evaluation of the performance of the power control algorithm proposed in Section IV. Note that the exact expressions of the SEs are used to compute the SEs. The simulation setup consists of a multi-cell network with 9 cells as illustrated in Fig. 2 addition, we use a wrap-around topology to avoid edge effects. We assume that the BSs are located in a 1 km 2 area. Each BS serves 5 CUs which are randomly distributed with uniform distribution in the BS's coverage area. The network contains 10 D2D pairs randomly distributed in the coverage area with uniform distribution. Since the pairs are distributed randomly, some cells will be more affected by D2D interference than others. We further assume that the distance between the transmitter and receiver of each D2D pair is 10 meters (note that the effect of D2D distance on spectral efficiency of CUs is further investigated in Appendix C). We use the three-slope pathloss model from [28] to model the large-scale fading in the network. This model contains two reference distances d 0 and d 1 , which are chosen to be 10 and 50 meters, respectively. The carrier frequency is 2 GHz, the bandwidth is 20 MHz [29], and the coherence interval contains 200 samples. The total number of pilots assigned for the D2D pairs is N = 5. First, we randomly select 5 D2D pairs and randomly assign a unique pilot sequence to each of them. For the remaining D2D transmitters, we randomly select and reuse a pilot sequence that is already used in the first group. Therefore, for the current setup, the minimum number of D2D pairs that is using a pilot sequence is one and the maximum is six. In addition, the noise variance is set to −94 dBm 3 and the maximum transmit power of the CUs and D2D transmitters are chosen to be P max = 200 mW for both pilot and data transmission.

A. Optimize data power control
In this subsection, we provide simulation results for data power control while all CUs and D2D transmitters use maximum power for pilot transmission. Figs. 3a and 3b show the cumulative distribution function (CDF) of the SE of a typical CU k for MR and ZF processing at the BSs, respectively, where the randomness is due to different user locations. We compare the proposed max-min data power control from Section IV-A and max product SINR power control from Section IV-A2 with three baseline methods: 4 • Equal power control (using full power); • Max-min data power control when there are only CUs; • Max product SINR power control when there are only CUs. It can be seen that the performance of our proposed maxmin and max product SINR algorithms are almost the same as in the corresponding case without underlaid D2D communication for both MR and ZF. Thus the power control can efficiently mitigate the D2D interference. For the maxmin power control, comparing the results with the case of maximum equal power transmission, max-min power control improves the performance of the 40% weakest users for MR and 0.064% of the weakest users for ZF processing. This is in line with the goal of max-min fairness power control, but it is clear that ZF by itself provides good performance for the most unfortunate CUs.
The proposed max product SINR power control algorithm improves the performance of 80% of users comparing with full power transmission when we have MR processing at the BSs. For the case of max product SINR with ZF processing at the BSs, there are only minor differences in SE in this case. Since ZF can efficiently mitigate the interference, making power control is less important in this case. This shows that with max product SINR, we get almost the same performance as with full power transmission. In fact, ZF mitigates the interference so well that everyone can transmit at full power in the network (Appendix D provides a further discussion on the effect of ZF processing on the data power coefficients of the D2D transmitters and CUs). From Fig. 3a and Fig. 3b it appears as if the SE of the system reduces when adding D2D users, but that is not the case. Fig. 4a and Fig. 4b, show the sum SE of all users in the network when using MR and ZF processing, respectively. It shows higher SE with D2D underlaying. This is due to the fact that we have more active users in the network with D2D underlay. Our proposed max product SINR power control shows a significant improvement in the sum SE for the case of D2D underlaying.

B. Optimized joint pilot and data power control
Next, we will jointly optimize the pilot and data powers. Fig. 5a and Fig. 5b show the SE of CU k for both cases of max-min and max product SINR power control with MR and ZF processing at the BSs, respectively. The figures compare the performance with joint pilot and data power control (denoted as "joint pilot & data PC") with the case of only data power control (denoted as "data PC"). Note that PC stands for power control in the figures. For both MR and ZF processing, the gain from joint pilot and data power control is larger    when having the max-min optimization objective. Note that the stopping criterion value ǫ ′ in Algorithm 2 is set to 0.001 in the simulations.
The sum SE is shown in Fig. 6a and Fig. 6b, for the cases of max-min and max product SINR power control with MR and ZF processing at the BSs, respectively. In both cases, the figures show that joint pilot and data power control for maxmin optimization objective enhances the SE of the network, as compared to only data power control. When we use max product SINR as the optimization objective, we have clear gains for MR processing which can be seen in Fig. 6a. But in the ZF processing case in Fig. 6b the curves are almost the same which reveals that ZF processing is capable of mitigating interference efficiently irrespectively of the power control. We observe from Fig. 5 and Fig. 6 that the gain from joint pilot and data power control, as compared to data power control only, is larger when using MR than with ZF. This is due to the fact that ZF can efficiently mitigate in interference spatially, if the channel estimates are good.

VI. CONCLUSION
In this paper, we have presented a framework for pilot allocation for multi-cell massive MIMO systems with underlaid D2D communication. Different orthogonal sets of pilots were used for cellular and D2D communication in order to mitigate extra interference at the cellular BSs. First, we derived the SE expressions for D2D pairs and the SE expression of CUs in the system for both MR or ZF processing at the BSs. Then, we proposed different power control schemes for the considered system model and solved max-min fairness and max product SINR optimization problems for MR and ZF processing at the BSs. We considered both data power control and joint pilot and data power control for both optimization objectives. Most problems were solved to global optimality, but the non-convex joint pilot and data power control with ZF processing at the BSs required an iterative convex approximation algorithm which only provides a local optimal point. The simulation results show that all the algorithms can effectively limit the interference between CUs and D2D pairs to enable D2D underlaying with limited performance degradation for the CUs. The joint pilot and data power control enhances the SE of the network in comparison with data power control for max-min optimization objective for both ZF and MR processing. Interestingly, the optimal power control with max product SINR optimization is rather similar to using full power. In addition, the proposed joint pilot and data power control provide larger SE gains in the case of MR processing comparing with ZF processing.

A. Tightness of the approximate SE for D2D communication
Since the SE expression for D2D communication was not in closed form, thus an approximation was provided in (35) to enable tractable power control optimization. Fig. 7 investigates the tightness of the approximation. We note that the approximation error is negligible for SEs less then 3 b/s/Hz, which is the range of main interest for max-min fairness optimization. For larger values of the SE, the approximation overestimates the SE by up to 0.5 b/s/Hz. Note that the D2D users are in relatively low SNR regime which is the reason for 0.5 b/s/Hz gap. When using the approximation for power control optimization, the D2D users will therefore get somewhat lower SE than what the output of the optimization problem predicts.

B. Basics of Geometric Programming
In this appendix, we provide the basics of geometric programming. We start by defining two types of functions. A monomial function is defined as where a > 0, the exponents b i are real-valued, and u is a vector where each element u i is the real and positive [30]. A posynomial function is formulated as the sum of one or multiple monomial functions and is written as where a k > 0 and the exponents b i,k are real-valued [30]. Adding or multiplying two posynomial functions produces a posynomial function as well. The result of dividing a posynomial function by a monomial function is also a posynomial. Dividing two posynomials result in a signomial function. This function has the same form as posynomial but the constants a k can be negative as well [30]. A geometric program on standard form is expressed as where the objective function f 0 (u) and the constraint functions f i (u) are posynomials, while the constraint functions g i (u) are monomial [30].

C. Effect of D2D distance on spectral efficiency
The distance between the transmitter and receiver of each D2D pair is one of the parameters which affects the SE of CUs. By increasing this distance, the interference from D2D transmitters to the CUs increases. The reason behind this is that D2D transmitters have to transmit with higher powers to compensate for the larger pathloss. Fig. 8 shows the effect of increasing the D2D distances on the SE of CU k with MR processing at the BSs. It shows that by increasing this distance from 10 m to 100 m, the SE of CU k decreases significantly which verifies the importance of the distance between the transmitter and receiver of D2D pairs on the SE of CUs. 10 m is in range of the typical distances between transmitter and receiver of D2D pairs which has been considered in the literature [12], [14], [18], [20]. The results in Fig. 8 verify that this distance is reasonable as it has a negligible effect on the SE of CU k. D. Effect of ZF processing on data power coefficients Fig. 9 shows the power coefficients of the D2D transmitters and CUs for max SINR product and max-min fairness with ZF processing at the BSs. The results show that for the case of max product SINR the D2D transmitters and CUs can transmit with full power as ZF can efficiently mitigate the interference. This is the reason behind the similarity between the performance of this case and full power transmission. Fig. 9: The power coefficients of the D2D transmitters and CUs for max SINR product and max-min fairness with ZF processing at the BSs.