QoS Aware Power Allocation and User Selection in Massive MIMO Underlay Cognitive Radio Networks

TL;DR

This paper presents a new optimization framework and a low-complexity algorithm for power allocation and user selection in massive MIMO underlay cognitive radio networks, aiming to maximize secondary users served while protecting primary users.

Contribution

It introduces the DMP algorithm for efficient user selection and power allocation, considering imperfect CSI and interference constraints, outperforming existing methods.

Findings

DMP algorithm serves more users with minimum rate constraints.

The approach approaches optimal performance with a large number of antennas.

Theoretical analysis quantifies interference and user rate exceedance.

Abstract

We address the problem of power allocation and secondary user (SU) selection in the downlink from a secondary base station (SBS) equipped with a large number of antennas in an underlay cognitive radio network. A new optimization framework is proposed in order to select the maximum number of SUs and compute power allocations in order to satisfy instantaneous rate or QoS requirements of SUs. The optimization framework also aims to restrict the interference to primary users (PUs) below a predefined threshold using available imperfect CSI at the SBS. In order to obtain a feasible solution for power allocation and user selection, we propose a low-complexity algorithm called DeleteSU-with-Maximum-Power-allocation (DMP). Theoretical analysis is provided to compute the interference to PUs and the number of SUs exceeding the required rate. The analysis and simulations show that the proposed DMP algorithm outperforms the stateof-the art selection algorithm in terms of serving more users with minimum rate constraints, and it approaches the optimal solution if the number of antennas is an order of magnitude greater than the number of users.

Figures

Equations

$${\mathbf{v}}_{\mathcal{S}(i)}^{\mathcal{S}}$$

$${\mathbf{G}}_{\mathcal{S}}$$

$$I_l=\sum _{k\in S}{I}_{kl}=\sum _{k\in S}{P}_k|{\mathbf{h}}_{l0}^H{\mathbf{v}}_k^{\mathcal{S}}{|}^2,l\in \mathcal{R},k\in \mathcal{S},$$

$$I_{jk}=P_j|{\mathbf{h}}_k^H{\mathbf{v}}_j^{\mathcal{S}}{|}^2,k,j\in \mathcal{S},j\ne k.$$

$$R_k^{\mathcal{S}}={\log }_2\left(1+\frac{P_k|{\mathbf{h}}_k^H{\mathbf{v}}_k^{\mathcal{S}}{|}^2}{{\sigma }_w^2+I_k+\sum _{j\in \mathcal{S},j\ne k}{I}_{jk}}\right),k\in \mathcal{S},$$

$${\hat{R}}_k^{\mathcal{S}}={\log }_2\left(1+\frac{P_k|{\hat{\mathbf{h}}}_k^H{\mathbf{v}}_k^{\mathcal{S}}{|}^2}{{\sigma }_w^2+I_k+{ϵ}_2}\right),k\in \mathcal{S}.$$

$$\underset{\{\mathcal{S},P_k\}}{\max }|\mathcal{S}|$$

$$:{\tilde{I}}_l=\sum _{k\in \mathcal{S}}{P}_k\left(|{\hat{\mathbf{h}}}_{l0}^H{\mathbf{v}}_k^{\mathcal{S}}{|}^2+{ϵ}_1\right)\le I^0,l\in \mathcal{R}$$

$${\hat{R}}_k^{\mathcal{S}}\ge R_k^0,k\in \mathcal{S},$$

$$\sum _{k\in \mathcal{S}}{P}_k\le P^0,P_k\ge 0.$$

$$\underset{\{\mathcal{S},P_k\}}{\max }|\mathcal{S}|$$

$$:\sum _{k\in \mathcal{S}}{P}_k{ϵ}_1\le I^0,$$

$$P_k$$

$$\sum _{k\in \mathcal{S}}{P}_k\le P^0.$$

$$(P1)\underset{\{s_k,P_k\}}{\max }$$

$$:\sum _{k=1}^K{s}_k{P}_k\le \min \left(I^0/{ϵ}_1,P^0\right),$$

$$P_k$$

$$s_k$$

$$P_k$$

$$:\sum _{k\in \mathcal{S}}{P}_k\le \min \left(I^0/{ϵ}_1,P^0\right),$$

$$P_k$$

$$P_k^{\mathcal{S}}=\frac{(2^{R_k^0}-1)\left({\sigma }_w^2+I_k+{ϵ}_2\right)}{|{\hat{\mathbf{h}}}_k^H{\mathbf{v}}_k^{\mathcal{S}}{|}^2},k\in \mathcal{S},$$

$$K_1^{\ast \ast }=\sum _{k\in {\mathcal{S}}_1^{\ast }}\text{mathbbm}1(R_k^{{\mathcal{S}}_1^{\ast }}\ge R_k^0),$$

$${\lambda }_k^{\mathcal{S}}=\frac{|{\hat{\mathbf{h}}}_k^H{\mathbf{v}}_k^{\mathcal{S}}{|}^2}{{\sigma }_w^2+I_k+{ϵ}_2},k\in \mathcal{S}.$$

$$P_k^{\mathcal{S}}=p_k^{\mathcal{S}}/{\lambda }_k^{\mathcal{S}},p_k^{\mathcal{S}}=(\mu {\lambda }_k^{\mathcal{S}}-1{)}^{+},k\in \mathcal{S},$$

$$\sum _{k\in \mathcal{S}}\left(\mu -\frac{1}{{\lambda }_k^{\mathcal{S}}}\right)=\min \left(\frac{I^0}{{ϵ}_1},P^0\right).$$

$$\hat{R}(\mathcal{S})=\sum _{k\in \mathcal{S}}{\log }_2\left(1+P_k^{\mathcal{S}}{\lambda }_k^{\mathcal{S}}\right).$$

$$\mathbb{E}[K_1^{\ast }]=\sum _{k=1}^{K}k\sum _{\mathcal{S}\in {\mathcal{S}}_k}f(\mathcal{S})g(\mathcal{S}),$$

$$f(\mathcal{S})$$

$$g(\mathcal{S})$$

$$P^{\prime }({\mathcal{S}}^{+}\backslash \{j\})$$

$$=(1-f({\mathcal{S}}^{+}))\underset{0}{\overset{\infty }{\int }}{p}_{P_j^{{\mathcal{S}}^{+}}}(x)\prod _{i\in {\mathcal{S}}^{+},i\ne j}{f}_{P_i^{{\mathcal{S}}^{+}}}(x)dx,$$

$${\kappa }_k^p$$

$${\theta }_k^p$$

$${\gamma }_k^{\mathcal{S}}=\frac{2^{R_k^0}-1}{({\beta }_k+{\sigma }_{\delta }^2)(M-|\mathcal{S}|-L+1)},$$

$${\gamma }_k^{{\mathcal{S}}_0}=\frac{2^{R_k^0}-1}{({\beta }_k+{\sigma }_{\delta }^2)(M-K-L+1)}.$$

$${\kappa }_p=\frac{{\left({\sum }_{j\in \mathcal{S}}{\kappa }_j^p{\gamma }_j^{\mathcal{S}}{\theta }_j^p\right)}^2}{{\sum }_{j\in \mathcal{S}}{\kappa }_j^p({\gamma }_j^{\mathcal{S}}{\theta }_j^p{)}^2},{\theta }_p=\frac{{\sum }_{j\in \mathcal{S}}{\kappa }_j^p({\gamma }_j^{\mathcal{S}}{\theta }_j^p{)}^2}{{\sum }_{j\in \mathcal{S}}{\kappa }_j^p{\gamma }_j^{\mathcal{S}}{\theta }_j^p},$$

$$\mathbb{E}[K_1^{\ast \ast }]=\sum _{k=1}^K\sum _{\mathcal{S}:k\in \mathcal{S}}f(\mathcal{S})g(\mathcal{S})\mathrm{Pr}(R_k^{\mathcal{S}}\ge R_k^0),$$

Full text

I Introduction

Due to increasing number of wireless devices and data rate demands, researchers are looking for various techniques to improve the spectrum efficiency of 5G wireless networks and serve a large number of devices in the available spectrum. Massive MIMO and underlay cognitive radio are being considered for 5G networks in order to accommodate more devices in the available spectrum [1, 2, 3]. In a massive MIMO system, a base station equipped with a large number of antennas serves multiple users using beamforming techniques in the same time-frequency resource block [4]. On the other hand, in an underlay cognitive radio (CR) network, a secondary base station (SBS) serves its users (secondary users) while keeping the interference to licensed primary users (PUs) below a specified threshold [5]. In underlay CR networks, the SBS transmits the downlink signal in the same time-frequency resource block as the primary transmitter. This is different from traditional interweave cognitive networks where the SBS transmits in an orthogonal time-frequency resource block.

The secondary BS, if equipped with a large number of antennas, can potentially employ beamforming techniques and serve multiple secondary users (SUs) in the downlink while limiting the interference to primary receivers. However, due to imperfect knowledge of the channels between PUs and the SBS, interference constraints at the PUs, and different rate requirements of SUs, it may not be feasible to serve all the SUs in the network [6]. Therefore, a judicious selection of SUs and power allocation are required at the SBS in order to simultaneously serve multiple SUs with required rates while limiting the interference to PUs.

I-A Related Work

Underlay CR networks with multiple antenna systems have received attention in recent years, since such networks allow concurrent transmissions from primary and secondary transmitters, thereby increasing the spectrum efficiency [7, 8, 9, 10, 11, 12, 13, 14, 15]. The works in [7, 8, 10, 11, 12, 14] consider small-scale MIMO with approximately ten or fewer antennas at the secondary transmitter. These works consider only one SU in the system, therefore they do not require SU selection. User selection mechanism has been partially considered in [9, 15]. An indirect selection mechanism is implemented in [9] where SUs receiving less than 0dB signal-to-interference-plus-noise-ratio (SINR) are dropped from the downlink transmission. The selection algorithm in [15] needs to know the number of users to be selected. A massive MIMO system has been employed in the secondary systems in [13, 16, 6]. In our previous works [16, 6], we proposed to use massive MIMO to serve multiple SUs concurrently with primary transmission, while the algorithm in [13] still serves only one SU. A selection algorithm under line-of-sight channels is proposed in [16], while the feasibility of serving all SUs under Rayleigh fading channels is studied in [6]. The selection algorithm was not considered in [6].

Massive MIMO systems differ from small-scale MIMO systems in [8, 9, 10, 14] in the design of beamforming (or precoding) vectors. In small-scale MIMO systems, optimum beamforming vectors are computed using iterative algorithms [8, 9, 14], and interference alignment [10]. Such approaches become prohibitively expensive in terms of complexity when used with massive MIMO systems. Using linear beamforming techniques such as zero-forcing (ZF), maximal-ratio combining (MRC) or minimum mean-square error (MMSE), the beamforming vectors can be computed using closed-form expressions without requiring any iterative search if the channels and the selected user set is known. In this paper, we focus on ZF beamforming as it is can also be used to restrict the interference toward PUs.

In an underlay CR network, the interference at primary receivers (PRs) can be eliminated using ZF beamforming if the channels between PRs and the SBS are perfectly known at the SBS. However, due to imperfect CSI in practical networks, there is non-zero leakage interference transmitted towards PRs even when ZF beamforming is used. The magnitude of the interference depends on the power allocated to SUs as well as the set of SUs selected. Therefore, there is a need to design a robust interference control mechanism along with power allocation and user selection in order to limit the interference to PRs below a specified threshold.

I-B Summary of Contributions and Outline

The main contributions of this work are summarized below.

A new optimization framework is proposed to select the maximum number of SUs in the downlink and obtain power allocation for the selected SUs in order to satisfy their instantaneous rate requirements. The interference to PRs is kept below a specified threshold using margin parameters that compensate for CSI estimation errors. The proposed formulation is different from the formulations in [17, 18, 19] which aim to maximize the sum-rate of the selected users and do not have interference constraints. 2. 2.

A new user selection and power allocation algorithm, called Delete-su-with-Maximum-Power-allocation (DMP), is proposed to obtain a feasible solution for the NP-hard optimization problem. Theoretical analysis of the algorithm is presented in order to compute the average number of SUs achieving the required rate, and average interference to primary receivers. The proposed algorithm is shown to achieve near-optimal results if the number of antennas at SBS is an order of magnitude larger than the number of users. 3. 3.

The user selection algorithm in [19] is extended for application in an underlay CR setting. The extended algorithm, called Modified Delete-Minimum-Lambda (MDML), also uses margin parameters and is robust against imperfect CSI. The proposed DMP selection algorithm is shown to serve more users than MDML in an underlay CR network.

Outline: This paper is organized as follows. The system model and the optimization problem are presented in Section II. The DMP and MDML algorithms are presented in Section III. Section IV presents the theoretical analysis and the optimality of the DMP algorithm. Simulation results are presented in Section V. Finally, the paper is concluded in Section VI.

Notations: We denote vectors by bold, lower-case letters, e.g., $\mathbf{h}$. Matrices are denoted by bold, upper case letters, e.g., $\mathbf{G}$. Scalars are denoted by non-bold letters e.g. $L$. Transpose, conjugate, and Hermitian of vectors and matrices are denoted by $(.)^{T}$, $(.)^{*}$, and $(.)^{H}$, respectively. The norm of a vector $\mathbf{h}$ is denoted by $||\mathbf{h}||$. $\Gamma(k,\theta)$ is the Gamma distribution with shape parameter $k$ and scale parameters $\theta$, whereas $\Gamma(x)$ is the Gamma function. The $i$-th element in the set $\mathcal{S}$ is denoted by $\mathcal{S}(i)$ and the cardinality of the set is denoted by $|\mathcal{S}|$. An empty set is denoted by $\emptyset$.

II System Model and Problem Formulation

II-A System Model

Consider an underlay CR network with one SBS and $K$ SUs. The SBS is equipped with $M(\gg K)$ antennas. This network coexists with $L$ primary transmitter-receiver pairs. Let $\mathcal{T}=\{1,2,...,L\}$ be the set of primary transmitters (PTs) and $\mathcal{R}=\{1,2,...,L\}$ be the set of PRs. The SUs, PRs and PTs are assumed to be single antenna terminals. Let $\mathbf{h}_{k}=\sqrt{\beta_{k}}\mathbf{\tilde{h}}_{k}\in\mathbb{C}^{M\times 1}$ be the channel between SU-$k$ and the SBS where $\beta_{k}$ is the slow fading coefficient accounting for attenuation and shadowing and $\mathbf{\tilde{h}}_{k}\sim\mathcal{CN}(0,\mathbf{I})$ [4, 20, 21]. The channel between PT-$l$ and SU-$k$ is denoted by ${h_{lk}}=\sqrt{\beta_{lk}}{\tilde{h}_{lk}}\in\mathbb{C},l\in\mathcal{T},k=\{1,2,...,K\}$, and ${\tilde{h}_{lk}}\sim\mathcal{CN}(0,1)$. Similarly, the channel between PR-$l$ and the SBS is $\mathbf{h}_{l0}=\sqrt{\beta_{l0}}\mathbf{\tilde{h}}_{l0}\in\mathbb{C}^{M\times 1},l\in\mathcal{R},\mathbf{\tilde{h}}_{l0}\sim\mathcal{CN}(0,\mathbf{I})$. We consider a time-division duplex (TDD) systems and the channels are assumed to be reciprocal. The network is depicted in Fig. 1.

The SBS has imperfect knowledge of the channels $\mathbf{h}_{l0},\mathbf{h}_{k}$. The estimates of channels are given by $\mathbf{\hat{h}}_{l0}=\mathbf{{h}}_{l0}+\mathbf{\Delta}_{l0}$ and $\mathbf{\hat{h}}_{k}=\mathbf{h}_{k}+\mathbf{\delta}_{k}$, respectively, where $\mathbf{\Delta}_{l0}\sim\mathcal{CN}(0,\sigma^{2}_{\Delta}\mathbf{I})$, and $\mathbf{\delta}_{k}\sim\mathcal{CN}(0,\sigma^{2}_{\delta}\mathbf{I})$ are the estimation errors. We model the quality of CSI between primary and secondary system using $\sigma^{2}_{\Delta}=\frac{\sigma^{2}_{w}}{P_{p}}$, while the quality of CSI within the secondary system is modeled using $\sigma^{2}_{\delta}=\frac{\sigma^{2}_{w}}{P^{0}}$, where $P_{p}$ and $P^{0}$ are transmit powers from PTs and the SBS, respectively, and $\sigma^{2}_{w}$ is the noise power at the SBS and SUs [22, 23, 24]. We consider block fading where the channels remain constant for a finite coherence interval.

Let $\mathcal{S}_{0}=\{1,2,...,K\}$ be the set of all SUs and $\mathcal{S}\subseteq\mathcal{S}_{0}$ be the set of SUs considered for downlink transmission and $P_{k},k\in\mathcal{S}$ be the power allocated to SU-$k$ in the downlink when set ${\mathcal{S}}$ is selected. The ZF transmit vector for SU-$k$ depends on $\mathcal{S}$ and is denoted by $\mathbf{v}^{\mathcal{S}}_{k}\in{\mathbb{C}^{M\times 1}},k\in\mathcal{S}$. The unit-norm ZF vectors are computed using [19]:

[TABLE]

where $\mathcal{S}(i)$ is the $i$-th entry in $\mathcal{S}$, and $[\mathbf{A}]_{i}$ is the $i$-th column of matrix $\mathbf{A}$. In the above equation, the matrix $\mathbf{G}_{\mathcal{S}}\in\mathbb{C}^{M\times(|\mathcal{S}|+L)}$ indicates the channel between SUs in set $\mathcal{S}$, PUs and the SBS. It should be noted that the ZF vector $\mathbf{v}^{\mathcal{S}}_{k}$ is in the null-space of the estimated channels to SU-$j,j\in\mathcal{S},j\neq k$. It is also in the null-space of the estimated channels to PRs. Therefore, $\mathbf{v}^{\mathcal{S}}_{k}$ satisfies $\mathbf{\hat{h}}^{H}_{j}\mathbf{v}^{\mathcal{S}}_{k}=0,j,k\in\mathcal{S},j\neq k$ and $\mathbf{\hat{h}}^{H}_{l0}\mathbf{v}^{\mathcal{S}}_{k}=0,l\in\mathcal{R}$.

The ZF beamforming vectors are not in the null-space of the true channel $\mathbf{h}^{H}_{l0}$ when the channel estimation error $\mathbf{\Delta}_{l0}$ is non-zero. Therefore, the interference caused at PR-$l,l\in\mathcal{R}$ is non-zero and can be expressed as:

[TABLE]

where $I_{kl}$ is the interference contribution of data stream of SU-$k$ towards PR-$l$. The interference $I_{l}$ depends on power allocation as well as the set of SUs selected. Similarly, the inter-SU interference at SU-$k$ due to the signal transmitted toward SU-$j$ can be expressed as

[TABLE]

Finally, the reverse interference at SU-$k$ is the sum of powers received from PTs: $I_{k}=\sum\limits_{l\in\mathcal{T}}P_{p}|{h_{lk}}|^{2}$, where $P_{p}$ is the power transmitted by PT-$l$.

II-B Optimization Problem

Our goal is to select the maximum number of SUs for downlink transmission in order to satisfy specific instantaneous rate of $R^{0}_{k}$ to selected SUs, while keeping the interference towards PRs below $I^{0}$. The total available power at the SBS is $P^{0}$. Note that the estimated interference to PR-$l,l\in\mathcal{R}$ based on the estimated channel is $\hat{I}_{l}=\sum_{k\in\mathcal{S}}P_{k}|\mathbf{\hat{h}}^{H}_{l0}\mathbf{v}^{\mathcal{S}}_{k}|^{2}=0$. Since the true interference $I_{l}$ is non-zero, we add a margin parameter $\epsilon_{1}$ to define $\tilde{I}_{l}=\sum_{k\in\mathcal{S}}P_{k}(|\mathbf{\hat{h}}^{H}_{l0}\mathbf{v}^{\mathcal{S}}_{k}|^{2}+\epsilon_{1})$ as the new estimated value of the interference with margin.

Further, the instantaneous rate achieved at SU-$k$, when a set $\mathcal{S}$ is selected, is

[TABLE]

where $\sigma^{2}_{w}$ is the noise power at the SU. Due to ZF beamforming, the estimated inter-SU interference will be zero, i.e., $\hat{I}_{jk}=P_{j}|\mathbf{\hat{h}}^{H}_{k}\mathbf{v}^{\mathcal{S}}_{j}|^{2}=0$ due to $\mathbf{\hat{h}}^{H}_{k}\mathbf{v}^{\mathcal{S}}_{j}=0,k\neq j$. We use a margin parameter $\epsilon_{2}$ to compensate for the estimation error in the channels between SBS and SUs. Therefore, the estimated instantaneous rate with margin becomes

[TABLE]

Note that, unlike $I_{l}$ and $I_{jk}$, the reverse interference term $I_{k}$ can be measured at SU-$k$ by observing combined signal received from all PUs during the channel estimation phase. The optimization problem can then be formulated as follows:

[TABLE]

Selection of parameters $\epsilon_{1}$ and $\epsilon_{2}$ is discussed in Section IV-B. By substituting $|\mathbf{\hat{h}}^{H}_{l0}\mathbf{v}^{\mathcal{S}}_{k}|^{2}=0$ and rearranging (8), we obtain the following equivalent optimization problem:

[TABLE]

Using binary selection variables $s_{k}\in\{0,1\}$ to indicate whether SU-$k$ is selected ($s_{k}=1$) or not ($s_{k}=0$), we can restate the above problem as:

[TABLE]

The constraint (15) is obtained by combining (11) and (13). This constraint indicates that the power allocation is controlled by the interference limit $I^{0}$ if $I^{0}/\epsilon_{1}\leq P^{0}$, while it is controlled by the transmit power limit $P^{0}$ if $I^{0}/\epsilon_{1}>P^{0}$.

III Selection algorithms and power allocation schemes

The optimization problem (14)-(17) to compute power allocations and selection variables is a non-convex mixed integer program and an NP-hard problem. Note that the computation of power allocations and selection variables depend on ZF vectors $\mathbf{v}^{\mathcal{S}}_{k}$ which in turn depend on the selected users. In order to solve this chicken-and-egg problem, we choose a particular set $\mathcal{S}$ and obtain ZF vectors and power allocation for that set. For a given selected set $\mathcal{S}$, the problem (14)-(17) reduces to the following feasibility problem with power allocations as variables:

[TABLE]

If the power allocation

[TABLE]

satisfies the constraint in (19), then it provides the solution to the power allocation problem for the set $\mathcal{S}$. Note that the above power allocation attempts to satisfy specific instantaneous rate of SU-$k$. Therefore, it is referred to as Qos-Aware-Power-allocation.

Let $K^{*}$ be the cardinality of optimal sets. The problem (14)-(17) can have multiple optimal sets, since multiple sets of the $K^{*}$ can satisfy the constraints (15)-(17). One approach of obtaining one of the optimal sets is to consider all possible sets of cardinalities $K,K-1,K-2,...,K^{*}$ one-by-one in decreasing order of cardinality, compute ZF vectors and power allocations by (1) and (21), respectively, and check whether the constraints in (19) are satisfied. Such approach of user selection is prohibitively complex and impractical since the number of sets to be considered increases exponentially with $K$. As an example, for $K=20$ and $K^{*}=5$, the minimum number of sets to be considered are $\sum_{K^{\prime}=K^{*}+1}^{K}$ ${K}\choose{K^{\prime}}$ $\approx 1$ million. Therefore, there is a need to design a low-complexity algorithm to select users and obtain power allocations.

As our goal is to maximize the cardinality of the set $\mathcal{S}$, we propose an approach which considers only one set of a particular cardinality that is obtained by dropping the SU that requires maximum power in a higher cardinality set. The selection algorithm is initialized by selecting all the SUs, i.e., $\mathcal{S}=\mathcal{S}_{0}$. ZF vectors $\mathbf{v}^{\mathcal{S}}_{k}$ and power allocations $P^{\mathcal{S}}_{k}$ are computed for the set $\mathcal{S}=\mathcal{S}_{0}$ using (1) and (21), respectively. Then, the condition $\sum_{k\in\mathcal{S}}P^{\mathcal{S}}_{k}\leq\min(I^{0}/\epsilon_{1},P^{0})$ is checked. If the condition is not satisfied, the SU with maximum power allocation is dropped from the set and a set $\mathcal{S}=\mathcal{S}\backslash\{j\}$ of lower cardinality is considered, where $j=\arg\max_{k\in\mathcal{S}}P^{\mathcal{S}}_{k}$. The ZF vectors and power allocations are re-computed for the new set using (1) and (21), respectively. This process is continued until the constraint $\sum_{k\in\mathcal{S}}P^{\mathcal{S}}_{k}\leq\min(I^{0}/\epsilon_{1},P^{0})$ is satisfied. Since the SU with maximum power allocation is dropped in each iteration, the algorithm is called Delete-su-with-Maximum-Power (DMP). Note that dropping of the SU that requires the maximum power causes maximum reduction in $\sum_{k}s_{k}P_{k}$ in constraint (15). This increases the probability that SUs included in set $\mathcal{S}\backslash\{j\}$ will satisfy the constraint (15). The algorithmic steps are summarized in Algorithm 1. It should be noted that the SUs which require excess power to satisfy their rate requirements will not be selected by the DMP. For example, if $P^{\mathcal{S}}_{k}=\frac{(2^{R^{0}_{k}}-1)\left(\sigma^{2}_{w}+{I}_{k}+\epsilon_{2}\right)}{|\mathbf{\hat{h}}^{H}_{k}\mathbf{v}^{\mathcal{S}}_{k}|^{2}}>\min(I^{0}/\epsilon_{2},P^{0})$, then SU-$k$ will not be selected.

The selected set is denoted by $\mathcal{S}_{1}^{*}$ and the cardinality of the selected set is $K^{*}_{1}=|\mathcal{S}^{*}_{1}|$. Due to imperfect CSI, all the selected SUs may not achieve the rate $R^{0}_{k}$. Therefore, we quantify the performance of the algorithm by $K^{**}_{1}(<K^{*}_{1})$ which is the number of SUs that achieve rate higher than $R^{0}_{k}$:

[TABLE]

where $\mathbbm{1}(.)$ is the indicator function.

We also consider a low-complexity version of this algorithm where step 9 of vector and power allocation update in Algorithm 1 is omitted. The set selected with this modified version without the vector is denoted by $\mathcal{S}_{2}^{*}$ and its cardinality by $K^{*}_{2}=|\mathcal{S}_{2}^{*}|$. Further, $K^{**}_{2}=\sum_{k\in\mathcal{S}^{*}_{2}}\mathbbm{1}(R^{\mathcal{S}^{*}_{2}}_{k}\geq R^{0}_{k}),$ is the number of SUs exceeding the required rate in DMP without vector update.

III-1 MDML Selection Algorithm

We extend the Delete-Minimum-Lambda (DML) selection scheme presented in [19] to underlay CR network with imperfect CSI. The DML algorithm in [19] selects users while maximizing the sum-rate of selected users. This approach does not take into account the rate constraints $R^{0}_{k}$ of SUs. It also does not take into account the imperfect CSI. Since DML was proposed for a primary massive MIMO network (CR network was not considered), it also does not include the reverse interference received at SU from primary transmitters. We modify the algorithm to include reverse interference and margin parameters for robustness against imperfect CSI. The Modified-DML (MDML) is described below.

In MDML, the equivalent channel gain between SU-$k,k\in\mathcal{S}$ and SBS is defined as:

[TABLE]

The power allocation for the SU included in the set $\mathcal{S}$ is obtained by water-filling. In order to satisfy the condition (15), the maximum power level of $\min\left({I^{0}/\epsilon_{1}},{P^{0}}\right)$ is used to compute the power allocation by water-filling as shown below

[TABLE]

where $(x)^{+}=\max\{x,0\}$, and $\mu$ is the water level satisfying

[TABLE]

Further, the estimated sum rate for SUs included in $\mathcal{S}$ can be written as

[TABLE]

The MDML algorithm drops the SU with minimum $\lambda^{\mathcal{S}}_{k}$, if dropping the SU results in increased sum-rate. The algorithmic steps are summarized in Algorithm 2. The selected set of SUs under this algorithm is denoted by $\mathcal{S}^{*}_{M}$.

IV Analysis of DMP algorithm

In this section, we provide analysis to compute $\mathbb{E}[K^{*}_{1}]$ and $\mathbb{E}[K^{**}_{1}]$ under fading channels $\mathbf{h}_{k},{h_{lk}}$ and $\mathbf{h}_{l0}$. We also analyze the average interference to primary receivers. Since the coefficients $\beta_{k},\beta_{lk},\beta_{l0}$ change slowly over time, they are assumed to be constant in the analysis [21].

IV-1 Average number of SUs served

The average number of SUs served using DMP is computed as follows:

[TABLE]

where $\mathcal{S}_{k}$ is the set of all sets of cardinality $k$, $f(\mathcal{S})$ is the probability of that the condition $\sum_{k\in S}P^{\mathcal{S}}_{k}\leq\min(I^{0}/\epsilon_{1},P^{0})$ is satisfied:

[TABLE]

and $g(\mathcal{S})$ is the probability of arriving at set $\mathcal{S}$ during the algorithmic iterations. Since the set $\mathcal{S}_{0}$ is always considered, we have $g(\mathcal{S}_{0})=1$ and $g(\mathcal{S}),\mathcal{S}\subset\mathcal{S}_{0}$ can be obtained using the following recursive expression:

[TABLE]

where $P^{\prime}\left(\{\mathcal{S}\cup\{j\}\}\backslash\{j\}\right)=P^{\prime}(\mathcal{S^{+}}\backslash\{j\})$ is the probability of dropping SU-$j$ from set $\mathcal{S^{+}}=\mathcal{S}\cup\{j\}$. This probability can be expressed as:

[TABLE]

where $p_{P^{\mathcal{S}^{+}}_{j}}(x)$ and $f_{P^{\mathcal{S}^{+}}_{i}}(x)$ are the probability density function (pdf) of $P^{\mathcal{S}^{+}}_{j}$ and the cumulative distribution function (cdf) of $P^{\mathcal{S}^{+}}_{i}$, respectively. In order to evaluate (28), we need distributions of $P^{\mathcal{S}}_{k}$ which can be obtained from Theorem 1. The distribution of $\sum_{k\in S}P^{\mathcal{S}_{k}}$ is required to evaluate (30) which can be obtained from Corollary 1.

Theorem 1**.**

The power allocation $P^{\mathcal{S}}_{k}$ in (21) is a Gamma random variable with shape and scale parameters $\kappa^{p}_{k}$ and $\theta^{p}_{k}$: $P^{\mathcal{S}}_{k}\sim\Gamma(\kappa^{p}_{k},\gamma^{\mathcal{S}}_{k}\theta^{p}_{k})$, where

[TABLE]

Similarly, for DMP without the vector update step, we have: $P^{\mathcal{S}_{0}}_{k}\sim\Gamma(\kappa^{p}_{k},\gamma^{\mathcal{S}_{0}}_{k}\theta^{p}_{k})$, where

[TABLE]

Proof.

Appendix A. ∎

Corollary 1**.**

Sum of powers $\sum_{k\in\mathcal{S}}P^{\mathcal{S}}_{k}$ follows the Gamma distribution: $\sum_{k\in S}P^{\mathcal{S}}_{k}\sim\Gamma(\kappa_{p},\theta_{p})$, where

[TABLE]

Proof.

Using Lemma 3 in [25], the sum $\sum_{k\in\mathcal{S}}P^{\mathcal{S}}_{k}$ is modeled as a Gamma random variable with shape and scale parameters $\kappa_{p}$, and $\theta_{p}$, respectively, as defined in (33). ∎

Corollary 2**.**

Consider selection of two sets $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ containing SU-$k$. The power required to achieve rate $R^{0}_{k}$ at SU-$k$ with selection of $\mathcal{S}_{1}$ stochastically dominates the power required to achieve the same rate with the selection of $\mathcal{S}_{2}$, if $\mathcal{S}_{2}$ is subset of $\mathcal{S}_{1}$, i.e., $\Pr(P^{\mathcal{S}_{1}}_{k}\geq x)>\Pr(P^{\mathcal{S}_{2}}_{k}\geq x)$ for any $x$, if $\mathcal{S}_{2}\subset\mathcal{S}_{1}$.

Proof.

Appendix B. ∎

Remark: In DMP, if set $\mathcal{S}_{1}$ does not satisfy constraints in (19), a subset of $\mathcal{S}_{1}$, say $\mathcal{S}_{2}$, is considered by dropping SU-$j$ that consumes the maximum power. Corollary 2 implies that the individual power requirements for all SUs still included in $\mathcal{S}_{2}$ reduce due to dropping of SU-$j$.

The expression for $\mathbb{E}[K^{*}_{2}]$ under DMP without vector update is obtained by evaluating $f(\mathcal{S})$ and $P^{\prime}({S^{+}\backslash\{j\}})$ using distributions of $P^{\mathcal{S}_{0}}_{k}$ instead of $P^{\mathcal{S}}_{k}$.

IV-2 Average number of SUs exceeding the required rate

The average number of SUs achieving the minimum rate of $R^{0}_{k}$ using DMP can be expressed as follows:

[TABLE]

In order to compute the above expression, we need to compute the complementary cdf of $R^{\mathcal{S}}_{k}$ which is obtained by Theorem 2.

Theorem 2**.**

The complementary cdf of the achieved rate at SU-k, if the set $\mathcal{S}$ is selected and the power is allocated by (21), is given by:

[TABLE]

where

[TABLE]

[TABLE]

Proof.

Appendix C. ∎

Remark: If $\sigma^{2}_{\delta}=0$, then we get $\Pr(R^{\mathcal{S}}_{k}\geq R^{0}_{k})=1$. The proof is provided in Appendix C.

The expression $\mathbb{E}[K^{**}_{2}]$ under DMP without vector update is obtained using the same expression as in the RHS of (34) by replacing $\gamma^{\mathcal{S}}_{j}$ with $\gamma^{\mathcal{S}_{0}}_{j}$ in Theorems 1 and 2.

IV-3 Average interference at PR-$l$

The expected value of the interference is computed as follows:

[TABLE]

where $\mathbb{E}\left[I_{kl}\right]=\mathbb{E}\left[P^{\mathcal{S}}_{k}|\mathbf{h}^{H}_{l0}\mathbf{v}^{\mathcal{S}}_{k}|^{2}\right]=\mathbb{E}\left[P^{\mathcal{S}}_{k}|\mathbf{\Delta}^{H}_{l0}\mathbf{v}^{\mathcal{S}}_{k}|^{2}\right]$. The second equality follows from $\mathbf{\hat{h}}^{H}_{l0}\mathbf{v}^{\mathcal{S}}_{k}=0$ due to zero forcing beamforming. The expression for $\mathbb{E}\left[I_{kl}\right]$ can be written using the distributions of $P^{\mathcal{S}}_{k}$ and $\mathbf{\Delta}_{l0}$ as follows:

[TABLE]

The proof is shown in Appendix D. Similarly, the expression for $\mathbb{E}[I_{l}]$ in DMP without vector update is obtained by substituting the following in (37):

[TABLE]

From (37), (38), and (39), we can see that the average interference to PRs is 0 for $\sigma^{2}_{\Delta}=0$.

IV-A Optimality of DMP

As described in Section III, the power allocated to SU-$k$, $P^{\mathcal{S}_{0}}_{k}$, remains constant during the algorithmic iterations of DMP without vector update. Therefore, DMP without vector update effectively obtains the solution for the following problem:

[TABLE]

The solution obtained by DMP without vector update can be written as follows:

[TABLE]

The solution $\mathcal{S}^{*}_{2}$ is the optimal solution for $\mathbf{P2}$, since no set of higher cardinality can satisfy the constraint (41) for fixed $P^{\mathcal{S}_{0}}_{k}$. This is because the proposed algorithm drops the SU with the highest power allocation in each iteration until the constraint (41) is satisfied and addition of any SU to the set $\mathcal{S}^{*}_{2}$ will violate the constraint.

Further, the optimal solutions of problems $\mathbf{P1}$ and $\mathbf{P2}$ differ due to the difference in power allocations $P^{\mathcal{S}^{*}}_{k}$ and $P^{\mathcal{S}_{0}}_{k}$, where $\mathcal{S}^{*}$ indicates the optimal solution set for $\mathbf{P1}$. The power allocations differ due to the difference in the number of nulls in the ZF vectors $\mathbf{v}^{\mathcal{S}^{*}}_{k}$ and $\mathbf{v}^{\mathcal{S}_{0}}_{k}$, which are denoted by $M-|\mathcal{S}^{*}|-L+1$ and $M-K-L+1$, respectively. The difference in the power allocations $P^{\mathcal{S}^{*}}_{k}$ and $P^{\mathcal{S}_{0}}_{k}$ becomes negligible if $M\gg K+L$. From Theorem 1, we can also observe that the distribution of $P^{\mathcal{S}^{*}}_{k}$ approaches that of $P^{\mathcal{S}_{0}}_{k}$ as $\gamma^{\mathcal{S}^{*}}_{k}\rightarrow\gamma^{\mathcal{S}_{0}}_{k}$, which occurs if $M\gg K+L$. Therefore, we can conclude that the problem $\mathbf{P1}$ becomes equivalent to $\mathbf{P2}$ for $M\gg K+L$ and $|\mathcal{S}^{*}_{2}|$ approaches the optimal solution $|\mathcal{S}^{*}|$.

Finally, the number of SUs selected by the DMP algorithm is no less than the number of SUs selected by DMP without vector update, i.e, $|\mathcal{S}_{2}^{*}|\leq|\mathcal{S}_{1}^{*}|\leq|\mathcal{S}^{*}|$. This is because the power allocated to each SU in DMP is less that or equal to that in DMP without vector update, i.e., $P^{\mathcal{S}}_{k}\leq P^{\mathcal{S}_{0}}_{k}$ when a set $\mathcal{S}$ is selected. This can also be observed from Corollary 2. The condition for selection of a set $\mathcal{S}$ under DMP $\left(\sum\limits_{k\in\mathcal{S}}P^{\mathcal{S}}_{k}\leq\min\left(\frac{I^{0}}{\epsilon_{1}},P^{0}\right)\right)$ is always satisfied if the condition under DMP without vector update is satisfied $\left(\sum\limits_{k\in\mathcal{S}}P^{\mathcal{S}_{0}}_{k}\leq\min\left(\frac{I^{0}}{\epsilon_{1}},P^{0}\right)\right)$, while the converse is not true. Therefore, we get $|\mathcal{S}_{2}^{*}|\leq|\mathcal{S}_{1}^{*}|\leq|\mathcal{S}^{*}|$. This phenomenon can be intuitively explained as follows. When the set $\mathcal{S}_{0}$ is selected initially, the number of degrees of freedom in the beamforming, after adding $L+K-1$ nulls, is $M-K-L+1$. The number of degrees of freedom increments if the ZF vectors are updated after dropping an SU. Therefore, the power requirements of SUs which are not dropped reduce as the ZF vectors $\mathbf{v}_{k}$ are better aligned with channels $\mathbf{\hat{h}}_{k}$. The reduction in power requirements implies that more SUs can be kept in the downlink, while satisfying the constraints of the problem $\mathbf{P1}$. Therefore, the number of SUs selected by the DMP is no less than the number of SUs selected by DMP without vector update.

IV-B Selection of algorithm parameters

The optimization framework $\mathbf{P1}$ in (14)-(17) involves various parameters. In this section, we provide discussion on the selection of these parameters. The parameters can be broadly classified into two categories: 1) network dependent fixed parameters: $P^{0},I^{0},R^{0}_{k}$, and 2) proposed margin parameters: $\epsilon_{1},\epsilon_{2}$.

IV-B1 Network dependent parameters

The network dependent parameters are decided by the operators of secondary and primary networks. Consider, for example, that primary and secondary networks coexist in 3.5GHz band as CBRS users where primary system is Priority Access License (PAL) user and secondary system is General Access Authorization (GAA) user [26, 27]. The value of $P^{0}$ will be determined using the power amplifier used at the SBS. Typical value of $P^{0}=40$dBm is used for BS under sub-6GHz bands. The rate constraints $R^{0}_{k}$ are determined by the operator of the secondary network depending on the QoS requirements for the SUs. The interference constraint $I^{0}$ is determined by the operator of primary network. The value of $I^{0}$ determines the SINR degradation of PUs due to the coexisting SUs. For example, if SINR degradation of $<1$dB is desired then $I^{0}$ should be set such that $I^{0}/\sigma^{2}_{w}<-6$dB, where $\sigma^{2}_{w}$ is the noise power at the PU.

IV-B2 Proposed margin parameters

Once the values of $P^{0},I^{0}$ and $R^{0}_{k}$ are fixed, the algorithm specific margin parameters $\epsilon_{1}$ and $\epsilon_{2}$ are set as described next. Margin parameters $\epsilon_{1}$ and $\epsilon_{2}$ are used in the proposed optimization framework in order to compensate for interference to PRs and inter-SU interference, respectively, resulting form imperfect CSI estimates.

In order to select appropriate value of $\epsilon_{1}$, let us consider average value of the true interference $I_{l}$ at PR-$l$ for given channel estimates $\mathbf{\hat{h}}_{l0}$ and a selected set $\mathcal{S}$:

[TABLE]

The last equality in the above equation follows from the fact that beamforming vectors $\mathbf{v}^{\mathcal{S}}_{k}$ are unit vectors that are in nullspace of estimated channels $\mathbf{\hat{h}}^{H}_{l0}$. Further, the proposed optimization problem ensures that $\sum_{k\in S}P^{\mathcal{S}}_{k}\epsilon_{1}\leq I^{0}$ due to constraint in (11). Therefore, the average interference $\mathbb{E}[I_{l}|\mathbf{\hat{h}}_{l0},\mathcal{S}]$ is below the threshold $I^{0}$ if $\epsilon_{1}\geq\sigma^{2}_{\Delta}$.

Now, let us consider the selection of $\epsilon_{2}$. From (4) and (5), we can see that the variable $\epsilon_{2}$ serves as a placeholder for inter-SU interference $\sum_{j\in\mathcal{S},j\neq k}I_{jk}$. The variable $\epsilon_{2}$ ensures that higher power is allocated to SU-$k$ to compensate for the inter-SU interference as seen from (21). Therefore, the value of $\epsilon_{2}$ should be selected such that $\epsilon_{2}\geq\sum_{j,j\neq k}I_{jk}$. However, the instantaneous value of inter-SU interference is unknown. We propose to set the parameter value such that it exceeds the expected value of inter-SU interference for given channel estimates $\mathbf{\hat{h}}_{j},\mathbf{\hat{h}}_{k}$ and a selected set $\mathcal{S}$:

[TABLE]

The equality $(a)$ results due to the fact that $\mathbf{\hat{h}}_{k}\mathbf{v}^{\mathcal{S}}_{j}=0$. The equality $(b)$ results from $\delta_{k}\sim CN(0,\sigma^{2}_{\delta}\mathbf{I})$ and $||\mathbf{v}^{\mathcal{S}}_{k}||=1$. Thus, $\epsilon_{2}$ should be greater than $\sum_{j\in{\mathcal{S}},j\neq k}P^{\mathcal{S}}_{j}\sigma^{2}_{\delta}$. However, the power allocations $P^{\mathcal{S}}_{j}$ are not known in advance. Therefore, we set $\epsilon_{2}\geq P^{0}\sigma^{2}_{\delta}$, which ensures that the condition in (46) is satisfied since $P^{0}\sigma^{2}_{\delta}\geq\sum_{j\in{\mathcal{S}},j\neq k}P^{\mathcal{S}}_{j}\sigma^{2}_{\delta}$.

IV-B3 Optimum margins $\epsilon_{1},\epsilon_{2}$

We observe that larger $\epsilon_{1}$ in (19) results in admitting fewer SUs in the downlink. Similarly, larger $\epsilon_{2}$ results in larger power allocation according to (21), further resulting in dropping of SUs due to the constraint $\sum_{k\in\mathcal{S}}P_{k}\leq\min(I^{0}/\epsilon_{1},P^{0})$ in (19). In order to satisfy the rate and interference constraints while admitting maximum number of SUs in the downlink, it is necessary to set $\epsilon_{1}$ and $\epsilon_{2}$ to the smallest possible values. Therefore, the setting $\epsilon_{1}=\sigma^{2}_{\Delta}$ and $\epsilon_{2}=P^{0}\sigma^{2}_{\delta}$ results in serving maximum number of SUs with given interference and rate constraints.

IV-C Complexity Analysis

Computational complexity of DMP as well as MDML is dominated by the computation of ZF vectors. For a set $\mathcal{S}$, the complexity of obtaining the ZF vectors is $\mathcal{O}\left(M(|\mathcal{S}|+L)^{3}\right)$ [19]. Since ZF vectors are updated in each iteration of DMP and MDML until a feasible set is reached, the worst case complexity of both algorithms is $\mathcal{O}\left(MK(K+L)^{3}\right)$, while the worst case complexity of DMP without vector update is $\mathcal{O}\left(M(K+L)^{3}\right)$.

As shown in the previous section, the solution $|\mathcal{S}^{*}_{2}|$ obtained by DMP without vector update approaches the optimal value with large $M$. Therefore, we can conclude that near-optimal number of SUs can be selected by the proposed algorithm while reducing the computational complexity by a factor of $K$ as compared to MDML.

V Simulation Results

In the results shown below, the noise power $\sigma^{2}_{w}$ is assumed to be $-100$dBm, the transmitted power from primary transmitters is $P_{p}=20$dBm, the transmit power limit is $P^{0}=40$dBm. The variance of error is modeled assuming reciprocal channels in a time-division duplex system as $\sigma^{2}_{\Delta}=\sigma^{2}_{w}/P_{p}$ and $\sigma^{2}_{\delta}=\sigma^{2}_{w}/P^{0}$ [22, 23, 24]. We consider uniformly distributed SUs and primary transmitters and receivers in a circular cell of radius $2$km with the SBS at the center. The minimum distance between the SBS and SUs is $100$m [28, 4]. For each realization of locations, the slow fading coefficients between two nodes are computed as $\beta=\rho d^{-3.8}$, where $d$ is the distance between the two nodes and $\rho$ is a log-normal shadowing variable with standard deviation $\sigma_{s}=8$dB. The margin parameters are set as $\epsilon_{1}=\sigma^{2}_{\Delta}$ and $\epsilon_{2}=P^{0}\sigma^{2}_{\delta}$. We simulate the algorithms for 1000 realizations of the channel for each realization of locations. Analytical and simulation results are averaged over 1000 realizations of the locations.

Comparison with optimal solution: The comparison of the average number of SUs selected by DMP and optimal selection is shown in Fig. 2. The optimal solution is obtained by considering all possible sets of cardinalities $K,K-1,K-2,...,K^{*}$ one-by-one in decreasing order of cardinality, computing ZF vectors and power allocations by (21), until the constraints in (19) are satisfied. We observe that the number of SUs selected by DMP is very similar to that by optimal selection and the difference between $E[|\mathcal{S}^{*}_{1}|]\approx E[|\mathcal{S}^{*}|]$. As the number of antennas increased from 64 to 256, the difference between the performance of DMP without vector update and optimal selection becomes negligible as explained in Section IV-A.

Impact of $R^{0}_{k}$: The proposed DMP algorithm is designed to satisfy the minimum rate required by the SUs unlike MDML which does not take into account the rate requirements. Therefore, it can be observed that the DMP serves more SUs exceeding the minimum required rate than MDML when $R^{0}_{k}$ is uniformly distributed in the range $(0,4]$ as seen in Fig. 3. The performance of the DMP and MDML becomes similar as the rate requirements increase to $4$bps/Hz. Further, it can be observed that the performance curves of the DMP and the MDML converge at a higher rate for large number of antennas. This indicates that the performance gain obtained by the DMP over MDML increases for a given rate requirement as the number of antennas increase.

Impact of $I^{0}$: The interference threshold $I^{0}$ limits the total transmitted power below $I^{0}/\epsilon_{1}$, thereby limiting the number of SUs served in both DMP and MDML. It should be noted that the interference of $-100$, $-106$ and $-110$ dBm results in SINR loss of $3,0.97$ and $0.41$dB, respectively at primary receivers. As shown in Fig. 4, the number of SUs served by the three algorithms increases by 1.5 times with increased interference threshold from $-110$dBm to $-100$dBm at the cost of reduced signal-to-interference-plus-noise ratio (SINR) at primary receivers. It can be observed that the performance gain obtained by the DMP over MDML is consistent for different interference thresholds.

Impact of number of primary tx-rx pairs $L$: Increased number of PTs in the network increases the reverse interference $I_{k}$ to SUs. This results in increased power requirement $P_{k}$ for the SU-$k$ to satisfy the rate according to (21). This increased power requirement in turn results in dropping of more SUs in step 6 of the DMP algorithm. Therefore, the number of SUs served by the proposed algorithm reduces with higher $L$ as shown in Fig. 5.

Impact of total number of SUs: The impact of increasing number SUs is shown Fig. 6. The number of SUs exceeding the required rate increases almost linearly for $M=128$ and $M=256$ under DMP when the rate requirements $R^{0}_{k}$ are uniformly distributed in the range $(0,4]$ bps/Hz. We also observe that the difference in the performance of the DMP with and without vector update reduces with increased number of antennas. This is due to the fact that the ratio $\gamma^{\mathcal{S}_{0}}/\gamma^{\mathcal{S}^{*}_{1}}$ is close to one which results in similar power allocations for SUs with and without vector update, thereby resulting in similar number of SUs being dropped in the step 6 of the DMP algorithm.

Optimality of margins: The margin parameters $\epsilon_{1}$ and $\epsilon_{2}$ are used to protect the PRs from the interference under imperfect CSI. In order to study the impact of margins on the performance, we plot the average interference and the number of SUs served for different values of $\epsilon_{1}$ and $\epsilon_{2}$. As shown in Fig. 7(a), the average interference remains below the threshold for $\epsilon_{1}>\sigma^{2}_{\Delta}$. This result holds for values of $\epsilon_{2}/P^{0}\sigma^{2}_{\delta}$ in range [0,4], because the variable $\epsilon_{2}$ does not significantly affect the average interference. As mentioned in Section IV-B, larger values of $\epsilon_{1}$ result in smaller the number of SUs served. Therefore, we keep $\epsilon_{1}=\sigma^{2}_{\Delta}=\sigma^{2}_{w}/P_{p}$ and plot $E[K^{**}]$ for range of values of $\epsilon_{2}$, as shown in Fig. 7(b). For $\epsilon_{2}<P^{0}\sigma^{2}_{\delta}$, fewer SUs receive the required rate due to inter-SU interference, while for $\epsilon_{2}>P^{0}\sigma^{2}_{\delta}$ fewer SUs are admitted in the downlink due to large power allocation. Therefore, we see that the maximum number of SUs are served for $\epsilon_{2}=P^{0}\sigma^{2}_{\delta}$ as described in Section IV-B3.

VI Conclusion

In this paper, we proposed an optimization framework in order to serve the maximum number of SUs in an underlay CR network consisting of a secondary BS equipped with a large number of antennas. The proposed framework uses margin parameters to limit the interference to PUs below a specified threshold under imperfect knowledge of CSI. A new user selection and power allocation algorithm, referred to as DMP, is proposed that is based on ZF beamforming and power allocation that satisfies specific rate requirements of selected SUs. Theoretical analysis is presented to compute the number of SUs selected and the interference caused at PUs by the proposed algorithm. Results show that the proposed DMP algorithm serves more SUs than modified DML algorithm for lower rate requirements. As the rate requirements for the SUs increase, the performance of the modified DML algorithm approaches that of DMP. A low complexity version of DMP without ZF vector update is also studied. This algorithm reduces the complexity by a factor of the number of SUs and provides similar performance as DMP with vector update when the number of SBS antennas is an order of magnitude larger than the number of SUs. The analysis and simulation results show that the number of SUs selected by the proposed algorithm approaches the optimal solution if the number of SBS antennas is an order of magnitude larger than the number of SUs and PUs in the network.

Appendix A Derivation of distributions of $P^{\mathcal{S}}_{k}$ and $P^{\mathcal{S}_{0}}_{k}$

From (21), $P^{\mathcal{S}}_{j}$ can be expressed as

[TABLE]

where $X=\sigma^{2}_{w}+{I}_{k}+\epsilon_{2}$ and $Y=|\mathbf{\hat{h}}^{H}_{k}\mathbf{v}^{\mathcal{S}}_{k}|^{2}=|({\mathbf{h}^{H}_{k}+\mathbf{\delta}^{H}_{k}})\mathbf{v}^{\mathcal{S}}_{k}|^{2}$. The vector $\mathbf{h}^{H}_{k}+\mathbf{\delta}^{H}_{k}\sim CN(0,\beta_{k}+\sigma^{2}_{\delta})$ is an isotropic vector, while the vector $\mathbf{v}^{\mathcal{S}}_{k}$ spans $M-|\mathcal{S}|-L+1$ dimensional space due to $|\mathcal{S}|-1+L$ nulls. Therefore, $Y$ is modeled as a Gamma random variable with shape and scale parameters $M-|\mathcal{S}|-L+1$ and $\beta_{k}+\sigma^{2}_{\delta}$, respectively [25, lemma 1] [6], i.e., $Y\sim\Gamma(M-|\mathcal{S}|-L+1,\beta_{k}+\sigma^{2}_{\delta})$. For simplicity of the analysis, we approximate $Y$ with its average value: $Y=(M-|\mathcal{S}|-L+1)(\beta_{k}+\sigma^{2}_{\delta})$. This approximation is valid because the variance of $Y$, $(M-|S|-L+1)(\beta_{k}+\sigma^{2}_{\delta})^{2}$, is negligible as compared to its mean $(M-|\mathcal{S}|-L+1)(\beta_{k}+\sigma^{2}_{\delta})$.

Further, consider the variable $X=\sigma^{2}_{w}+{I}_{k}+\epsilon_{2}$. We define a constant $C=\sigma^{2}_{w}+\epsilon_{2}$. Each term in the summation ${I}_{k}=\sum_{l\in\mathcal{T}}P_{p}|{{h}_{lk}}|^{2}$ is modeled as a Gamma random variable with distribution $\Gamma(1,P_{p}\beta_{lk})$. Therefore, the mean of $X$ is $C+\sum_{l\in\mathcal{T}}P_{p}\beta_{lk}$, while its variance is $\sum_{l\in\mathcal{T}}{P_{p}}^{2}\beta_{lk}^{2}$. We model $X$ as a Gamma random variable with size and shape parameters $\kappa^{p}_{k}$ and $\theta^{p}_{k}$, respectively, i.e., $X\sim\Gamma(\kappa^{p}_{k},\theta^{p}_{k})$ Therefore, we have

[TABLE]

[TABLE]

By solving for $\kappa^{p}_{k}$ and $\theta^{p}_{k}$, we obtain (31).

Finally, $P^{\mathcal{S}}_{k}=\frac{2^{R^{0}_{k}}-1}{(M-|S|-L+1)(\beta_{k}+\sigma^{2}_{\delta})}X=\gamma^{\mathcal{S}}_{k}X$ is a Gamma random variable with size and shape parameters $\kappa^{p}_{k}$ and $\gamma^{\mathcal{S}}_{k}\theta^{p}_{k}$, respectively, i.e., $P^{\mathcal{S}}_{k}\sim\Gamma(\kappa^{p}_{k},\gamma^{\mathcal{S}}_{k}\theta^{p}_{k})$. Note that $P^{\mathcal{S}}_{k},k=1,2,...$ are independent variables since they are functions of independent random variables ${{h}_{lk}}$.

The distribution of $P^{\mathcal{S}_{0}}_{k}$ is obtained by following the above derivation with $Y=(M-K-L+1)(\beta_{k}+\sigma^{2}_{\delta})$. This is due to the fact that the vector $\mathbf{v}^{\mathcal{S}_{0}}_{k}$, in this case, spans $M-K-L+1$ dimensional space due to $K-1+L$ nulls.

Appendix B Proof of corollary 2

Consider two sets $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ containing SU-$k$ such that $\mathcal{S}_{2}\subset\mathcal{S}_{1}$ and $|\mathcal{S}_{1}|>|\mathcal{S}_{2}|$. Since $P^{\mathcal{S}_{1}}_{k}$ and $P^{\mathcal{S}_{2}}_{k}$ are Gamma random variables, their CDFs can be written as follows:

[TABLE]

Since $\gamma^{\mathcal{S}_{1}}_{k}>\gamma^{\mathcal{S}_{2}}_{k}$, we get $\Pr(P^{\mathcal{S}_{1}}_{k}\leq x)<\Pr(P^{\mathcal{S}_{2}}_{k}\leq x)$ or $\Pr(P^{\mathcal{S}_{1}}_{k}\geq x)>\Pr(P^{\mathcal{S}_{2}}_{k}\geq x)$.

Appendix C Derivation of $\Pr(R^{\mathcal{S}}_{k}\geq y)$

The rate achieved at SU-$k$ is given as

[TABLE]

where $\mathcal{S}$ is the selected set. Substituting for $P_{k}^{\mathcal{S}}$ from (21) in the above equation, we get

[TABLE]

Similar to the variable $Y$ in the previous appendix, variables $|\mathbf{h}_{k}^{H}\mathbf{v}^{\mathcal{S}}_{k}|^{2}$ and $|\mathbf{\hat{h}}_{k}^{H}\mathbf{v}^{\mathcal{S}}_{k}|^{2}$ are approximated with their average values $(M-|\mathcal{S}|-L+1)\beta_{k}$ and $(M-|\mathcal{S}|-L+1)(\beta_{k}+\sigma^{2}_{\delta})$, respectively. Let us define $C_{y}$ as follows:

[TABLE]

Substituting the above equation, ${I}_{k}=\sum_{l}P_{p}|{{h}_{lk}}|^{2}$, and $I_{jk}=P_{j}|\mathbf{h}_{k}^{H}\mathbf{v}_{j}|^{2}$ in (51), we can rewrite the equation as follows:

[TABLE]

where $Z_{lk}=(1-C_{y})P_{p}|{h_{lk}}|^{2},Z_{k}=\sum_{j\in\mathcal{S},j\neq k}P^{\mathcal{S}}_{j}|\mathbf{h}_{k}^{H}\mathbf{v}^{\mathcal{S}}_{j}|^{2}$, and $\zeta_{y}$ as defined in (36). Since $Z_{lk}$ and $Z_{k}$ are independent random variables, the cdf in the RHS of (53) can be expressed in terms of Fourier transforms of the characteristic functions these variables. Therefore, we derive the characteristic functions of $Z_{lk}$ and $Z_{k}$. The variable $Z_{lk}$ is a Gamma random variable $\sim\Gamma(1,(1-C_{y})P_{p}\beta_{lk})$ with characteristic function:

[TABLE]

where $\theta^{z}_{lk}=(1-C_{y})P_{p}\beta_{lk}$. Further, the variable $Z_{k}=\sum_{j\in\mathcal{S},j\neq k}P^{\mathcal{S}}_{j}|\mathbf{h}_{k}^{H}\mathbf{v}^{\mathcal{S}}_{j}|^{2}$ can be written as $Z_{k}=\sum\limits_{j\in\mathcal{S},j\neq k}P^{\mathcal{S}}|\mathbf{\hat{h}}_{k}^{H}\mathbf{v}^{\mathcal{S}}_{j}+\mathbf{\delta}_{k}^{H}\mathbf{v}^{\mathcal{S}}_{j}|^{2}$ $=\sum_{j\in\mathcal{S},j\neq k}P^{\mathcal{S}}_{j}|\mathbf{\delta}_{k}^{H}\mathbf{v}^{\mathcal{S}}_{j}|^{2}$. The second equality follows from $\mathbf{\hat{h}_{k}v^{\mathcal{S}}_{j}}=0$ due to zero forcing beamforming. The term $|\mathbf{\delta}_{k}^{H}\mathbf{v}^{\mathcal{S}}_{j}|^{2}$ is the projection of isotropic vector $\delta_{k}\sim CN(0,\sigma^{2}_{\delta})$ on uncorrelated space spanned by $\mathbf{v}^{\mathcal{S}}_{k}$, which gives $|\mathbf{\delta}_{k}^{H}\mathbf{v}^{\mathcal{S}}_{j}|^{2}\sim\Gamma(1,\sigma^{2}_{\delta})$ [25, lemma 3]. Therefore, $P^{\mathcal{S}}_{j}|\mathbf{\delta}_{k}^{H}\mathbf{v}^{\mathcal{S}}_{j}|^{2}$ is a product of two Gamma random variables and is approximated as a Gamma random variable [29]. The mean and the variance of $P^{\mathcal{S}}_{j}|\mathbf{\delta}_{k}^{H}\mathbf{v}^{\mathcal{S}}_{j}|^{2}$ are given below:

[TABLE]

The variable $Z_{k}$ is modeled as a Gamma random variable with shape parameter $\kappa^{z}_{k}$ and shape scale parameter $\theta^{z}_{k}$. The parameters are computed using moment matching method [25, lemma 3] by solving the following two equations for $\kappa^{z}_{k}$ and $\theta^{z}_{k}$:

[TABLE]

Expressions in (36) follow from (55) and (56). The characteristic function of the Gamma random variable $Z_{k}$ is as follows [30]:

[TABLE]

Since $Z_{lk}$ and $Z_{k}$ are independent random variables, the cdf in the RHS of (53) can be written in terms of the Fourier transform of the product of characteristic functions of these random variables as follows:

[TABLE]

Substituting for $\phi_{lk}(jt)$ and $\phi_{k}(jt)$ from (54) and (57), respectively, we obtain the expression (35). The expression for $\Pr(R^{\mathcal{S}}_{k}\geq y)$ under DMP without vector update is computed by following the above derivation and replacing $\gamma^{\mathcal{S}}_{j}$ by $\gamma^{\mathcal{S}_{0}}_{j}$ in (55).

If we have $\sigma^{2}_{\delta}=0$, $\mathbf{\hat{h}}_{k}=\mathbf{{h}}_{k}$, and $I_{jk}=0$. Then, substituting $y=R^{0}_{k}$ in (51), we get $\Pr(R^{\mathcal{S}}_{k}\geq R^{0}_{k})=1$.

Appendix D Derivation of $\mathbb{E}[I_{lk}]$

In DMP, we have $\mathbb{E}\left[I_{kl}\right]=\mathbb{E}[P^{\mathcal{S}}_{k}|\mathbf{\Delta}^{H}_{l0}\mathbf{v}^{\mathcal{S}}_{k}|^{2}]$. The expression $\mathbb{E}[P^{\mathcal{S}}_{k}|\mathbf{\Delta}^{H}_{l0}\mathbf{v}^{\mathcal{S}}_{k}|^{2}]$ is obtained by following the derivation of $\mathbb{E}[P^{\mathcal{S}}_{j}|\mathbf{\delta}_{k}^{H}\mathbf{v}^{\mathcal{S}}_{j}|^{2}]$ in Appendix C and replacing $\mathbf{\delta}_{k}$ and $\sigma^{2}_{\delta}$ with $\mathbf{\Delta}_{l0}$ and $\sigma^{2}_{\Delta}$, respectively. In DMP without vector update, the expression for $\mathbb{E}\left[I_{kl}\right]$ is obtained by replacing $\gamma^{\mathcal{S}}_{k}$ with $\gamma^{\mathcal{S}_{0}}_{k}$ in (38).