Worst-Case Energy Efficiency Maximization in a 5G Massive MIMO-NOMA System

In this paper, we examine the robust beamforming design to tackle the energy efficiency (EE) maximization problem in a 5G massive multiple-input multiple-output (MIMO)-non-orthogonal multiple access (NOMA) downlink system with imperfect channel state information (CSI) at the base station. A novel joint user pairing and dynamic power allocation (JUPDPA) algorithm is proposed to minimize the inter user interference and also to enhance the fairness between the users. This work assumes imperfect CSI by adding uncertainties to channel matrices with worst-case model, i.e., ellipsoidal uncertainty model (EUM). A fractional non-convex optimization problem is formulated to maximize the EE subject to the transmit power constraints and the minimum rate requirement for the cell edge user. The designed problem is difficult to solve due to its nonlinear fractional objective function. We firstly employ the properties of fractional programming to transform the non-convex problem into its equivalent parametric form. Then, an efficient iterative algorithm is proposed established on the constrained concave-convex procedure (CCCP) that solves and achieves convergence to a stationary point of the above problem. Finally, Dinkelbach’s algorithm is employed to determine the maximum energy efficiency. Comprehensive numerical results illustrate that the proposed scheme attains higher worst-case energy efficiency as compared with the existing NOMA schemes and the conventional orthogonal multiple access (OMA) scheme.

significant chunk of network cost. EE is defined as the amount of information transmitted per joule of energy consumed [13,14]. In addition, combining NOMA together with MIMO can further increase the SE and EE performance in future wireless networks [15][16][17][18][19]. A significant number of researchers has examined the performance of the NOMA scheme with the distinct topics such as user pairing algorithms, power allocation strategies, beamforming and fairness schemes in order to study the feasibility of adopting the NOMA scheme for 5G communications. The beamforming (BF) technique can minimize the inter beam interference (IBI) by allocating a single beamforming vector to a group of users. Thus, the integration of NOMA with multi-user BF (NOMA-BF) can further enhance the system capacity in future wireless systems [20,21].
However, in order to reduce the inter user interference (IUI), the impact of user pairing was examined for a MIMO-NOMA system in [22,23], where the authors exploit the fixed power allocation (FPA) scheme and pair the users with distinctive channel conditions assuming perfect channel state information (CSI) at the base station (BS). Proportional fairness (PF) based user pairing algorithm is proposed in [24,25] to maximize the achievable rate for the single cell downlink NOMA system. A greedy user pairing and iterative power allocation (PA) algorithm was proposed in [26] to maximize the throughput in an orthogonal frequency division multiplexing (OFDM) based NOMA system. A hierarchical PA scheme is proposed in [27] based on the distinctive channel gains to maximize the sum rate and expedite a bigger number of users in a NOMA system. The authors in [28] proposed a distributed matching algorithm to optimize the user pairing and the transmit power allocation between the near and far users. The proposed user pairing algorithm helps to minimize the inter-user interference (IUI) and maximize the system throughput in a downlink NOMA network. A novel user pairing scheme is proposed in [29] for 5G cellular network, where the users having the largest proportional fairness (PF) coefficient is paired with the users having the best channel conditions in order to mitigate the IUI. A dynamic user scheduling and clustering design is proposed in [30] based on the limited feedback information at the BS to suppress the IUI and inter-cluster interference in the downlink space division multiple access NOMA system. Moreover, the proposed algorithm reduces the system overhead and maximizes the net throughput. The physical layer security of a single-input single-output (SISO) NOMA system was investigated by authors in [31], where they proposed a optimal power allocation strategy to maximize the secrecy sum rate (SSR) and successive interference cancellation (SIC) technique is employed at the receiver side to cancel the IUI. In [32], the authors proposed a low complexity user pairing scheme based on the greedy search method to maximize the weighted sum-rate (WSR) in a downlink NOMA system. A virtual user pairing scheme is proposed in [33], where the multiple near users having a similar channel gain is paired with a single far user to effectively utilize the spectrum in the NOMA system.
We examine the EE maximization problem for multicell massive MIMO-NOMA downlink system with imperfect channel state information at the transmitter (CSIT). Here, we especially consider an ellipsoidal uncertainty model (EUM) [27] to form the CSI errors, where the robust design of the beamforming matrix is taken into account to solve the worst-case EE maximization problem. Since NOMA is an interference bound scheme, we also propose a novel joint user pairing and dynamic power allocation (JUPDPA) algorithm, which mitigates the inter user interference (IUI) and also helps to attain the maximum EE. Several approaches for the EE maximization problem have been proposed in the literature [34][35][36][37] for the perfect CSI at the BSs. For instance, authors in [34] proposed an energy efficient optimal PA technique for a SISO-NOMA system that maximizes the EE subject to the minimum data rate requirement for each user. EE resource allocation problem has been studied in [35], where they proposed a sub optimal resource allocation algorithm to optimize the sub channel and PA factor in order to maximize the EE for a single cell NOMA system. In [36], the authors investigated the EE-SE co-design for the single cell and multi-cell NOMA system considering the perfect CSI at the BS. Successive convex approximation (SCA) based approach is proposed in [37] to tackle the EE maximization problem for a cognitive radio inspired multiuser downlink MISO-NOMA system, which maximizes the EE with respect to the quality of service (QOS) constraint for each primary user.
However, in practice, it is generally hard to achieve the perfect CSI in fading channels. Thus, it is meaningful to deal with the imperfect CSI in the problem formulation.
Few related works in [38,39] studied the EE maximization problem considering the imperfect CSI at the transmitter. In particular, authors in [38] examined the EE maximization problem for a MIMO-NOMA system with statistical CSI at the transmitter side and proposed a near optimal PA scheme to maximize the EE. An EE resource allocation algorithm was proposed in [39] to maximize the secret EE and secret key EE for both MISO-single antenna eavesdropper and MIMO-multiple antennas eavesdropper systems considering perfect and statistical CSI at the base station. Although there exists an enormous interest in EE maximization problem, most of the authors considered the problem with the assumption of perfect CSIT or statistical CSIT. Thus, it largely remains unsolved in typical scenarios of interest. For instance, a worst-case beamforming design to maximize the EE for multicell massive MIMO-NOMA system with imperfect CSIT has not been examined before. Motivated from the above observations, the major contribution of this work is summarized as follows.

•
In this paper, we investigate the EE maximization problem in a multicell massive MIMO-NOMA downlink system. The EE maximization problem is formulated under the constraints of transmission power and the minimum signal-to-interference-plus-noise-power-ratio (SINR) requirement for the CE users. The considered problem has non-convex fractional objective function, which is intractable and difficult to solve.

•
We first turn the fractional objective function into an equivalent parametric one by introducing a non-negative parameter. Successive convex approximation (SCA) method is further employed to solve the dual problem for minimizing the transmission power at the BS and also to find the upper bound values for the original convex optimization problem.

•
We propose a joint robust beamforming and EE maximization (JRBEE) algorithm established on the constrained concave-convex procedure (CCCP) and the Dinkelbach method that solve and achieve the convergence to a stationary point of the above problem. The JRBEE algorithm determines the optimal beamforming vector, which helps to maximize the worst-case EE. • Moreover, we also propose a novel JUPDPA algorithm based on the median and the Euclidean norm of the channel vectors. The JUPDPA algorithm helps to minimize the inter-user interference and also enhance the fairness between the users by dynamically allocating the transmission power to the paired users.

•
Simulation results guarantee the convergence property of our proposed EE maximization algorithm and showed that the proposed scheme can significantly improve the worst-case energy efficiency compared to the existing NOMA schemes and the conventional OMA scheme.

Organization
The remainder of this paper is organized as follows. Section 2 describes the considered system model and formulates the worst-case EE maximization problem. A novel JUPDPA scheme is presented in Section 3. A robust beamforming algorithm using CCCP is proposed to tackle the EE maximization problem and find the optimal EE in Section 4. Comprehensive numerical results are presented to show the excellent performance of our proposed scheme in Section 5, and Section 6 concludes this paper.

System Model and Problem Formulation
We examine the downlink transmission in a multicell massive MIMO-NOMA system as shown in  Multicell Massive multiple-input multiple-output non-orthogonal multiple access (MIMO-NOMA) downlink system.

Transmitter
The transmitted signal vector x from n B antennas is given as x = VF, where V = [V 1,n , V 2,n , V 3,n , .., V M,n ] and F = F 1,n , F 2,n , F 3,n , .., F M,n are the beamforming matrix and the encoded signal for all the M pairs, respectively. The BS transmits the symbol vector T 1,m,n and T 2,m,n in the same time and frequency slot via power domain NOMA. The superimposed signal of the cell center (CC) and cell edge (CE) users belonging to the mth pair in nth cell is written as F m,n = P 1,m,n T 1,m,n + P 2,m,n T 2,m,n , where T 1,m,n and T 2,m,n are the desired signal for the CC and CE user, respectively. P 1,m,n and P 2,m,n are the PA coefficients for the CC and CE user, respectively, which have to satisfy P 1,m,n + P 2,m,n = 1 for each user pair. V m,n ∈ C 1×n M,m represents the beamforming vector corresponding to pair m, where m ∈ M. The transmit power constraints for each BS is given by where P p is the maximum power transmitted at the BS and p ∈ N B .

Channel Model
The overall channel matrix for all the user pairs can be expressed as follows: where H ∈ C n M ×n B and M is the total number of user pairs in a single cell. The channel matrix corresponding to m th user pair in n th cell is written as where Ψ m,n in Equation (4) represents the CSIT uncertainty, which is considered to lie in the bounded ellipsoidal region as where m ∈ M, n ∈ N and G 0 is a given matrix which represents the shape of the region and ε 2 manages the size of the ellipsoidal region. In this paper, we consider that G is of full rank so that H m,n has a structural definition of being an ellipsoid [27].Ĥ m,n is the minimum mean square error (MMSE) estimate of channel matrix which is obtained from [40] aŝ where p m,n is the average transmit power of each uplink pilot symbol and B m,n is a Gaussian matrix with i.i.d. CN (0, 1) entries. The users present in each cell simultaneously transmit training sequences of length φ m,n symbols at the start of every coherence interval.

Receiver
The received signal vector at the mth pair in nth cell is given as y m,n = y 1,m,n + y 2,m,n , where y 1,m,n and y 2,m,n are the received signal by CC and CE users, respectively.

Cell Center User
The received signal by the CC user is given as y 1,m,n = h T 1,m,n V m,n P 1,m,n T 1,m,n + M ∑ j=1,j =m h T 1,m,n V j,n F j,n + z 1,m,n , where F j,n = P 1,j,n T 1,j,n + P 2,j,n T 2,j,n and z 1,m,n is the additive white gaussian noise (AWGN) noise with CN (0, σ 2 1,m,n ). In Equation (8), the first term represents the desired signal to a CC user and the second term is due to the interference from other pairs present in the same cell, which we call the inter pair interference (IPI). Inter user interference (IUI) that occurs due to CE users present in the same pair is cancelled due to the implementation of SIC at the CC user. The SINR corresponding to the CC user is represented as

Cell Edge User
The received signal by the CE user is given as where Ξ = h T 2,m,n V m,n P 1,m,n T 1,m,n , Π = ∑ M j=1,j =m h T 2,m,n V j,n F j,n and z 2,m,n is the AWGN noise vector with CN (0, σ 2 2,m,n ). Variable g T 2,m,c represents the interfering channel vector (from other cell) to the CE user present at the mth pair of nth cell. In Equation (10), the first term is the desired signal to CE user and the second term is the interference from CC user present in the same pair which we call it as IUI. The third term is the IPI and the fourth term is due to the interference from the CE users present in other cells, which is called inter cell interference (ICI). The SINR corresponding to the CE user is represented as

Energy Efficiency
The achievable downlink rate of the transmission from BS to m th pair is the mutual information between the unknown transmitted signal F m,n and the observed received signal y m,n , and is obtained as where R 1,m,n and R 2,m,n are the achievable downlink rates by the CC and CE users, respectively. The achievable rates for CC and CE users are given in Equations (13) and (14), respectively, as The sum-rate for all pairs is expressed as where R 1,m,n and R 2,m,n as given in Equations (13) and (14), respectively. The total transmit power at the base station is given as where P p is the maximum transmission power at the BS and P c is the circuit power consumption at each BS. Therefore, the energy efficiency (EE) can be defined as the ratio of overall sum rate to the total power consumption, which is expressed as

Problem Formulation
The energy efficiency maximization problem for massive MIMO-NOMA system with ellipsoidal uncertainty model (EUM) is formulated as follows where EE is defined in Equation (17) and Γ thd 2,m,n is the SINR threshold for the CE users. The objective function in (18) maximizes the energy efficiency with bounded channel uncertainties. The constraints given in (18) represent the minimum SINR requirement for the CE users and the transmit power constraint for each BS.

User Pairing Scheme
In this section, we propose a novel user pairing scheme that reduces the inter-user interference (IUI) between the paired users and also enhance the fairness between them. Algorithm 1 summarizes the steps of the proposed joint user pairing and dynamic power allocation (JUPDPA) algorithm. In Step 1, users are ordered in the descending order based on their Euclidean norm of the channel vector ĥ k,m,n and then we split the available users into two categories, i.e., CC and CE with the help of the median given in (19) for the effective user pairing. We define the estimated channel vector from here on as asĥ x instead ofĥ T x,m,n for simplicity. The median (MED) for the total available users in a cell can be expressed as After splitting, we extricate the users present in both CC and CE group as the odd and even numbered users. The proposed user pairing scheme is given in Step 2, where the odd users from the CC group select the least difference Euclidean norm user from the CE group and even users from the CC group choose the CE user, which has the large difference in the Euclidean norm between them. In Step 3, we eliminate the paired users in Step 2 from the CC and CE user sets and repeat the above steps until CC or CE user set is empty. The transmission power for the paired users will be allocated based on the proposed dynamic power allocation (DPA) scheme in Step 4, while the unpaired users obtain their desired signal power with the help of the conventional OMA scheme in Step 5. Algorithm halts in the final step.

Dynamic Power Allocation Scheme
Once the user pairing is done, we need to allocate the PA coefficients to the CC and CE user. PA is done dynamically for each selected user pair in order to further maximize the EE of the MIMO-NOMA system within each user pair. Before doing that, the PA limit for the CC user needs to be fixed and the difference between the PA for each CC user has to be calculated, which is given as The PA coefficient for the user closest to the BS in the CC group starts from L (derived in (24)) and it will allocate 2L to the next nearest user. This process goes on until the maximum PA limit for the CC user is reached. Once we set the PA coefficient for all the CC users, the power allotted to CE user can be found from P 2,m,n = 1 − P 1,m,n . The main advantage of the proposed JUPDPA scheme is that the fairness between the CC and CE users are considerably improved due to the increase in the achievable rate of the CE users. In particular, the achievable rate of CE users paired with odd CC user will get the huge benefits from our proposed scheme. Although the achievable rate of few CE users paired with even CC users are slightly degraded due to their least difference Euclidean norm pairing, but it is acceptable considering the improvement in the fairness between the users and also the overall system performance. This enhancement is mainly due to the maximum PA limit for the CC users. In general, energy-efficient transmission strategy transmits data at a very low rate in order to reduce the power consumption. SINR for the CC users is relatively high due to their presence close to the base station (BS) and therefore less transmission power is sufficient to satisfy their data rate requirements. Thus, in our problem, we consider the SINR constraint only on the CE users, which are present far way from the BS and also suffer from low SINR due to the IUI as well as the inter-cell interference (ICI) from the neighboring cells. The proposed JUPDPA algorithm pairs a CC user with the CE user to reduce the IUI. The proposed JUPDPA algorithm also allocates more than 2/3 of the overall transmission power at the BS to CE users in order to meet their data rate requirements and also to improve the system fairness. Moreover, the proposed JRBEE maximization algorithm in Section 4 also helps to attain the data rate requirement for the CE users by finding the optimal beamformers for each pair (CC and CE user), which, in turn, minimize the inter pair interference (IPI) and the power consumption at the BS. The transmit power constraint given in the original problem (18) for each BS helps to mitigate the interference from the neighboring cells to the CE users. It is good to note that the worst-case EE is set to zero, if the maximum transmission power is not large enough to meet the SINR threshold value for the CE users.
Step 1. Initialization and splitting: where D represents the Euclidean norm of channel vectors for all the available users in the descending order and k ∈ N M . The users in D are categorized into CC and CE group as given in (20) and (21), respectively, based on the above condition and the median (MED) given in (19): where i and j represents the CC and CE user index, respectively, ( i ∈ I and j ∈ J ).
Step 2. User Pairing: Select a user from CC group and calculate the difference in the Euclidean norm between that user and all the available users from the CE group, which are ordered in descending order. The proposed user pairing scheme is given in the below figure.
If the selected user from the CC group is odd numbered, then choose the least difference Euclidean norm user from the CE group for pairing. However, for the even numbered CC user, choose the user from the CE group that has the largest difference in the Euclidean norm of the channel vectors.

Algorithm 1 Cont.
Step 3. Remove and Repeat: Remove the selected users H m,n = [ ĥ i , ĥ j ] from (20) and (21). Save them in H = [H 1,n , H 2,n , .., H M,n ], where m represents the user pair index: Step 4. Power Allocation: Paired users will be allocated transmission power based on the dynamic power allocation scheme as given in Table 1.
Step 5. Unpaired Users: Remaining unpaired users will be serviced based on the OMA scheme.
Step 6. Exit: Halt the Algorithm. Table 1 represents the dynamic PA scheme assuming five CC users, while it is good to note that the proposed scheme can be generalized for any number of CC users following a similar procedure. Based on the above observations, the power allotted to the nearest CC user and farthest CE will be 0.06 and 0.94, respectively, considering five users in the CC group. The overall PA for all the CC users (considered in Algorithm 1) are given in Table 1.

Proposed Joint Robust Beamforming and Energy Efficiency Maximization Algorithm
Here, in this section, we propose an efficient iterative algorithm based on the CCCP [41,42] and the Dinkelbach algorithm to find the optimal beamforming vector (V * m,n ) for each paired user and determine the maximum energy efficiency, respectively. The Dinkelbach algorithm converts the fractional non-convex optimization problem into a sequence of auxiliary problems while ensuring a super linear convergence rate with limited complexity [43,44]. The worst-case EE maximization problem formulated in (18) is intractable and non-deterministic polynomial hard (NP-hard). Moreover, it is very difficult to solve due to its nonlinear fractional objective function and semi infinite non-convex SINR constraint of the CE users given in (18). In order to solve the above NP-hard problem, we first implement the Dinkelbach algorithm to convert the fractional objective function into a parametric equivalent one with the help of the non-negative parameter λ. The original problem given in (18) is reformulated as follows: where The objective function of the reformulated problem is still nonlinear, since it is neither concave nor convex and it comes under the category of the difference of convex problem. We follow the CCCP procedure to solve the difference of convex programming problem given in (25). CCCP mainly helps to convert the non-convex problem into a sequence of convex problems by iteratively finding a linear approximation of the original non-convex objective function around the current solution and then solve the approximated convex function in each iteration. In addition, CCCP also guarantees convergence to the stationary point of the original convex optimization problem [41,42]. In the process of converting the above problem into a tractable one, we first define the transmit covariance matrices L m,n V m,n V † m,n 0, since the objective function and constraints present in the original problem (18) are linear in the matrix VV † . The optimization problem after substituting the change of variables is as follows where EE (λ, L m,n ) = R sum ({L m,n } M m=1 ) − λ P T and the achievable rate for CC user is written as where f 1 (L m,n ) is a concave function and g 1 (L m,n ) is a convex function, which are given below in (29) and (30), respectively, as The achievable rate for the CE user is given as where f 2 (L m,n ) is a concave function and g 2 (L m,n ) is a convex function, which are given below in (32) and (33), respectively, as where A ∑ c∈N \n ∑ M m=1 g T 2,m,c V m,c 2 . Its good to observe that we have employed the semi definite relaxation (SDR) method to relax the rank constraint, i.e., rank{L m,n } = 1 with the semi definite matrix L m,n 0 [45]. The non-convex problem in (27) is well approximated by the convex objective function given below in (35). First order Taylor approximation helps to replace g 1 (L m,n ) and g 2 (L m,n ) with their convex majorants by the following inequality that supports due to the concavity of log(x) as Z(U, C) log det(C) + 1 ln tr(C −1 (U − C)).
Considering the above changes, the optimization problem in (27) can be reformulated as follows: where Θ = ∑ M j=1,j =m h T 2,m,n V j,n 2 and ∆ = ∑ c∈N \n ∑ M m=1ĝ 2,m,c (L m,c )ĝ † 2,m,c . The problem in (35) is still non-convex due to the presence of infinitely many SINR constraints occurred due to the channel uncertainties. Therefore, we need to convert the semi infinite non-convex SINR constraint of the CE user into the linear convex constraint. Fortunately, in literature, we have a well known S-procedure that gives a linear matrix inequality (LMI) representation for adaptive quadratic constraint over an uncertainty set characterized by quadratic inequalities. Equivalent LMI representation for the constraints in (35) are given by converting those infinite constraints into a finite number of linear constraints with the help of auxiliary variables α m for m ∈ M. For convenience, we cite the S-procedure in [46] here. Let P and Q be two n × n Hermitian matrices, r ∈ C n and s ∈ R. Then, the two following conditions are equivalent: (a) a H Pa + r H a + a H r + s ≥ 0 for all a ∈ C n such that a H Qa ≤ 1; (b) there exists a d ∈ R such that d ≥ 0 and By straightforward implementation of S-procedure, the constraint given in (35) is transformed as I n B ⊗ Ω m,n + α m,n I n B ⊗ G m,n vec(ĥ 2,m,n Ω m,n ) vec(Ω m,nĥ2,m,n ) † χ 0, where α m,n ≥ 0, m = {1, 2, 3, .., M}, χ = tr(ĥ 2,m,n Ω m,nĥ † 2,m,n ) − Γ thd 2,m,n I − α m,n ε 2 , Ω m,n = P 2,m,n L m,n − Γ thd 2,m,n (P 1,m,n L m,n + where B = ∑ c∈N \n ∑ M m=1ĝ 2,m,c (L m,c )ĝ † 2,m,c . It is worth mentioning here that the uncertainty present in the transformed constraint is now dealt by the parameters i.e., G m,n and ε. Now, problem (35) can be equivalently reformulated to a convex optimization problem as given below: where M m ({L m,n } M m=1 ) = I n B ⊗ Ω m,n + α m,n I n B ⊗ G m,n vec(ĥ m,n Ω m,n ) vec(Ω m,nĥm,n ) † tr(ĥ 2,m,n Ω m,nĥ † 2,m,n ) − Γ thd 2,m,n I − α m,n ε 2 . (41) Now, we observe that (40) is a tractable convex problem that can be solved with the help of convex packages available in the literature [47]. However, we can find the sub optimal beamforming vector for each pair via solving the problem in (40), and the total transmit power is increased to maintain the required SINR for each CE user. Thus, in order to enhance the accuracy of first order Taylor approximation and to find the upper bound values for the original convex optimization problem, and we consider the dual problem of power minimization at the base station under the SINR constraint for the CE user. The dual optimization problem is formulated subject to the SINR constraints for all the CE users is as follows: where T 1,m,n =ĥ † 1,m,n L m,nĥ1,m,n , T 2,m,n =ĥ † 2,m,n L m,nĥ2,m,n and The considered dual problem is non-convex due to the infinite number of SINR constraints on the CE user. To solve the above problem, we employ the successive convex approximation (SCA) approach [47] and further implement S-Lemma to convert the infinite number of SINR constraints into a finite one. The sub-optimal beamforming vectors obtained by solving the convex problem in (40) is considered as the initial values for the transmit power minimization problem. The convex reformulated problem for minimizing the transmission power is given as where M m ({L m,n } M m=1 ) are already given in (41). Algorithm 2 summarizes the steps of the proposed joint robust beamforming and energy efficiency (JRBEE) maximization algorithm where the optimal beamforming vectors are obtained iteratively to maximize the EE. It is good to note that the non-decreasing optimal values attained over the proposed JRBEE maximization algorithm can minimize the power consumption at the BS. Since both (45) and (46) are convex, updating L m,n iteratively will only increase or maintain the objective value, which converges to a stationary solution of the original problem (35). The optimal beamforming vectors V * m,n found with the proposed JRBEE maximization algorithm are substituted in (17) to maximize the EE by finding the optimal value i.e., λ * , where it yields EE(λ * , V * m,n ) = 0. In the proposed JRBEE maximization algorithm, we first initialize {L (k) m,n } M m=1 to arbitrary positive semidefinite matrix, which is feasible for the problem (27). The initialization value plays a key role in deciding the convergence rate of the proposed JRBEE maximization algorithm. If the chosen arbitrary initialization value is very close to the optimal solution, then it takes only fewer iterations to attain the stationary point. Otherwise, if the chosen arbitrary initial value is not close to the optimal solution, then the proposed iterative algorithm will still achieve the same objective value as the previous one, but the number of iterations taken for the convergence will increase. We then initialize the matrix {L (l) m,n } M m=1 for the transmit power minimization problem (46) by using the beamforming vector (V m ) obtained from solving the problem (45) in the proposed JRBEE maximization algorithm. We also initialize λ init = 4 for the outer loop iteration, which employs Dinkelbach method to determine the maximum EE. Moreover, Figure 3 in Section 5 proves that the proposed energy efficient maximization scheme converge to a stationary point within few iterations using the arbitrary feasible initialization point.
until convergence criterion is satisfied. 3: end for 2.2 Calculate the linear beamformer (V m ) for m pairs and SINR for all users.
until convergence criterion is satisfied. 3: end for 2.4 Determine the optimal beamformer (V * m,n ) that achieves minimum transmission power.

3.
Obtain the EE by substituting V * m,n into the (26) and find the value of λ. 4. Update λ in step 1 and repeat the procedure until the convergence criteria is satisfied to find λ * .

Complexity Analysis for the Proposed Algorithms
We first calculate the complexity of the proposed JUPDPA algorithm, where all the available users are split into two categories, i.e., CC and CE users. Each pair contains only two users to avoid the SIC complexity at the receiver side. CC users chose their corresponding paired user based on the proposed JUPDPA algorithm by exhaustively searching all the available users in the CE category.
Thus, the computational complexity of the proposed JUPDPA algorithm is O( and J m is the total number of available users and the number of users present in each pair (we assume two users per pair), respectively. The computation complexity of the proposed JRBEE maximization algorithm mainly consists of two parts. The first part is the required number of iterations taken for the convergence and the second part is the complexity taken to solve the convex optimization problem in (45) and (46) during each iteration. The overall computational complexity for the proposed JRBEE maximization algorithm is O(Z 3 ((Z 1 + Z 2 )(N M N B ) 3.5 logµ −1 )), where Z 1 and Z 2 are the two inner loop iterations required to solve the convex problem (45) and (46), respectively. Z 3 is the number of outer loop iterations required to determine the maximum EE using the Dinkelbach method, which is known for its super linear convergence rate. N B and µ is the number of BSs and the solution accuracy, respectively. The convex (CVX) tool is used to solve the above convex optimization problems using interior point approach, which takes O((N M N B ) 3.5 logµ −1 ) computations [48]. The proposed energy efficiency maximization (EEmax) scheme (JUPDPA + JRBEE maximization algorithm) has the least implementation complexity due to the simpler matrix calculations. Average run-time computed in Figure 4 also proves that the proposed EEmax scheme takes less time to attain the objective value. The beamforming matrix is determined as V m,n ← B m,n C 1/2 m,n , where the diagonal elements present in C m,n matrix represents the non-zero eigenvalues of L D m,n (D ∈ {k, l}) and the corresponding eigenvectors are present in the column space of the matrix B m,n .

Numerical Results and Discussion
In this section, we evaluate the performance of our proposed EEmax scheme (includes JUPDPA and JRBEE algorithm) and compare it with the existing NOMA schemes in [20,21] and the conventional OMA scheme [49]. The proposed JUPDPA in Algorithm 1 and JRBEE algorithm given in Algorithm 2 are jointly labeled as proposed EEmax Scheme. The compared existing schemes are as follows: (1) we obtain the EE considering EUM when the sum rate is maximized for the iterative beamforming algorithm proposed in [20] referred to as Ratemax scheme in the sequel; (2) user pairing and PA algorithm along with the zero-forcing beamforming (ZFBF) proposed in [21] is evaluated for EUM to achieve the EE referred to as EE-ZFBF scheme. Moreover, the OMA transmission with imperfect CSI is considered where two users transmit in two orthogonal time slots and the robust beamforming design is utilized to tackle the corresponding optimization problem via CCCP to maximize the EE, labeled as OMA scheme [49]. We consider a multi cell massive MIMO-NOMA downlink system as shown in Figure 1. We assume there are three cells, N B = 3 base stations, each with multiple transmit antennas and N M users, each equipped with a single receive antenna. Each BS is subject to the same transmit power constraint. Unless otherwise specified, we assume channel uncertainty (error bound) for EUM is ε = 0.05, β = 0 dB, Γ thd 2,m,n = 10 dB, N M = 10, n B = 20, λ init = 4, P P = 30 dBm, P C = 1 W, and σ 2 m,n = 1. We will examine the worst-case energy efficiency (EE) performance for EUM. Worst-case EE was achieved by averaging the energy efficiency over the acquisition of fading channel matrices. Unless explicitly stated, all results are presented for 20 channel fading realizations. Jain's fairness index [50] is employed to measure the fairness between the users. We also presume that the elements of the channel matrices between the user in the b th cell and BS in the a th cell follows independent and identically distributed (i.i.d.) complex Gaussian with CN (0, β |a−b| ) distribution in which we call β the inter-cell channel gain. We consider only two users in each pair to reduce the SIC receiver's decoding complexity. Convergence criteria is set to 10 −4 , i.e., we assume that the proposed JRBEE algorithm comes to halt when the difference between attained worst-case EE values within the successive iterations are less than 10 −4 . Figure 2 demonstrates the worst-case energy efficiency versus available transmission power at the BS with the fixed circuit power of P C = 1 W. It can be observed that the proposed EEmax scheme consistently outperforms the Ratemax and OMA scheme in terms of the worst-case EE. It can be noticed that the worst-case EE of the proposed scheme increases until the transmission power reaches 28 dBm and begins to saturate after that level. This can be explained considering that the BS stops consuming more power to transmit the signals when the proposed EEmax scheme achieves the optimal EE for the given convergence criteria. In particular at P P = 30 dBm, the worst-case EE achieved by the proposed EEmax scheme is 17.6 Mbits/joule compared to 16.5 Mbits/joule by the EE-ZFBF scheme. The worst-case EE for the Ratemax scheme increases until the power reaches 28 dBm and then it begins to decrease, since the scheme emphasis on finding the BF vectors to achieve the high transmission capacity. It is also good to note that the worst-case EE performance of the proposed EEmax scheme for EUM performs close to the case with perfect CSI.  Figure 3 illustrates the worst-case EE versus number of iterations for two different transmission powers at the BS. It can be noticed that the worst-case EE increases monotonically for the proposed EEmax scheme before it converges at around five iterations for both low (P P = 15 dBm) and high (P P = 30 dBm) transmission power. However, the Ratemax and OMA scheme takes around eight and ten iterations, respectively, for convergence. This is due to the reason that the optimal BF vectors in the proposed JRBEE algorithm are obtained by minimizing the transmission power at the base station. We also note that the proposed EEmax scheme consistently beats the Ratemax scheme in terms of the worst-case EE and the number of iterations to attain convergence. Figure 4 depicts the number of antennas at the transmitter versus an average CPU run time, which helps to determine the implementation complexity of the considered schemes. It can be seen that the proposed EEmax scheme has the least implementation complexity compared to the Ratemax and the OMA scheme. In particular, the average runtime taken by the proposed scheme for 20 antennas is only 65 s, whereas the Ratemax scheme and the OMA scheme takes around 90 s and 135 s, respectively. This observation can be justified for the reason that the proposed user pairing algorithm splits the available N M users into two groups and search only N M 2 times to allocate the transmission power for selecting one active user pair. Moreover, the proposed EEmax scheme also avoids the large matrix inversions and complex matrix calculations.     Figure 5 investigates the worst-case EE versus SINR requirement for the CE user. It can be observed that the proposed EEmax scheme consistently outperforms the EE-ZFBF and Ratemax schemes in terms of the worst-case EE. The worst-case EE gradually decreases for all the schemes while there is a increment in the minimum SINR requirement for the CE users. In particular, at SINR = 20 dB, the worst-case EE attained by the proposed EEmax scheme is 13.4 Mbits/joule compared to the 12.5 Mbits/joule and 10.5 Mbits/joule achieved by the EE-ZFBF and the Ratemax scheme, respectively. This performance gain can be validated in considering that the proposed EEmax scheme allocates more than 2/3 of the total available transmission power to the CE users. Moreover, the BF vectors attained by the proposed JRBEE maximization algorithm are robust even in the worst-case conditions. It is also good to note that the worst-case EE is set to zero if the transmit power is not large enough to meet the SINR requirement of the CE users. Figure 6 exemplifies the average fairness index versus transmission power at the BS for a different number of users. We compare our proposed JUPDPA algorithm with RUPFPA (random user pairing with FPA) [22] and PFUPPA (proportional fairness based user pairing and PA) [23] algorithm. It can be noticed that the fairness linearly increases with the transmission power with the proposed JUPDPA algorithm providing the best fairness among the considered schemes. In particular, while considering for 10 users at P P = 25 dBm, the average fairness index is 0.83 for the proposed JUDPA algorithm compared to 0.79 and 0.7 for the PFUPPA and RUPFPA algorithm, respectively. This was expected because the proposed pairing scheme has a significantly greater impact on the achievable rate of the CE users, which, in turn, helps to increase the fairness between the users. Moreover, it is good to note that the average fairness index is reduced for a large number of users due to increase in the IUI.   In Figure 7, the worst-case energy efficiency is evaluated versus the number of users for different channel uncertainty values. It can be noticed that when the number of users increases, the worst-case energy efficiency also increases only until all the users needs can be satisfied by the total available power at the BS. Once it exceeds that limit, then the worst-case EE will start to decrease. Worst-case EE attained by the proposed EEmax scheme decreases with an increase of channel uncertainty as expected. Our proposed EEmax scheme again achieves better performance than the other considered schemes. In particular, at = 0.09 for 10 users, the worst-case energy efficiency for the proposed EEmax scheme is 11.5 Mbits/joule compared to 8.5 Mbits/joule attained by the OMA scheme. This can be justified by the fact that the conventional OMA scheme does not fully utilize the available spectrum and energy resources compared to NOMA scheme where a greater number of users can be served simultaneously by the power domain multiplexing.

Conclusions
In this paper, a robust beamforming design was studied to solve the EE maximization problem in a 5G massive MIMO-NOMA system with imperfect CSI at the BS. In particular, we have considered the EUM for the inclusion of CSI errors. In order to tackle the above problem, we first converted the fractional nonlinear objective function into its equivalent parametric form. An efficient iterative algorithm based on the CCCP was proposed to optimize the transmit beamforming vectors for each pair. We then employed the Dinkelbach method to find the optimal parameter that maximize the worst-case EE. We also proposed a novel JUPDPA algorithm that improves the fairness between the users and maximize the worst-case EE by dynamically allocating the transmission power to the paired users. Numerical results guarantee that the proposed EEmax scheme converges faster than the Ratemax scheme and the OMA scheme. Via numerical results, it was also confirmed that the worst-case EE achieved by the proposed EEmax scheme persistently outperforms the existing NOMA schemes and the conventional OMA scheme.

Abbreviations
The following abbreviations are used in this manuscript: