1

Resource Optimization with Load Coupling
in Multi-cell NOMA
Lei You∗ , Di Yuan∗ , Lei Lei† , Sumei Sun‡ , Symeon Chatzinotas† , and Björn Ottersten†

arXiv:1708.05281v3 [] 7 Jan 2018

∗ Department of Information Technology, Uppsala University, Sweden
{lei.you; di.yuan}@it.uu.se
† Interdisciplinary Centre for Security, Reliability and Trust, Luxembourg University, Luxembourg
{lei.lei; symeon.chatzinotas; bjorn.ottersten}@uni.lu
‡ Institute for Infocomm Research, A*STAR, Singapore
sunsm@i2r.a-star.edu.sg

Abstract—Optimizing non-orthogonal multiple access (NOMA)
in multi-cell scenarios is much more challenging than the singlecell case because inter-cell interference must be considered. Most
papers addressing NOMA consider a single cell. We take a
significant step of analyzing NOMA in multi-cell scenarios. We
explore the potential of NOMA networks in achieving optimal
resource utilization with arbitrary topologies. Towards this goal,
we investigate a broad class of problems consisting in optimizing
power allocation and user pairing for any cost function that is
monotonically increasing in time-frequency resource consumption. We propose an algorithm that achieves global optimality for
this problem class. The basic idea is to prove that solving the joint
optimization problem of power allocation, user pair selection,
and time-frequency resource allocation amounts to solving a
so-called iterated function without a closed form. We prove
that the algorithm approaches optimality with fast convergence.
Numerically, we evaluate and demonstrate the performance of
NOMA for multi-cell scenarios in terms of resource efficiency
and load balancing.
Index Terms—NOMA, multi-cell, resource allocation

I. I NTRODUCTION

N

ON-ORTHOGONAL multiple access (NOMA) is considered as a promising technique for enhancing resource
efficiency [2]–[13]. In two recent surveys [2], [3], the authors
pointed out that resource allocation in multi-cell NOMA
poses much more research challenges compared to the singlecell case, because optimizing NOMA with multiple cells
has to model the interplay between successive interference
cancellation (SIC) and inter-cell interference. As one step
forward, the investigations in [2], [3] have addressed two-cell
scenarios. In [6], the authors proposed two coordinated NOMA
beamforming methods for two-cell scenarios. Reference [8]
uses stochastic geometry to model the inter-cell interference in
NOMA. Hence the results do not apply for analyzing network
with specific given network topology. Reference [9] optimizes
energy efficiency in multi-cell NOMA with downlink power
control. However, the aspect of determining which users share
resource by SIC, i.e., user pairing, is not considered. To the
best of our knowledge, finding optimal power allocation and
user pairing simultaneously for enhancing network resource
efficiency in multi-cell NOMA without restrictions on network
topology has not been addressed yet.
Part of this paper has been presented at IEEE GLOBECOM, Singapore,
Dec. 2017 [1].

The crucial aspect of multi-cell NOMA consists of capturing the mutual interference among cells; This is a key
consideration in SIC of NOMA. Therefore, the cells cannot
be optimized independently. For orthogonal multiple access
(OMA) networks, a modeling approach had been proposed
that characterizes the inter-cell interference via capturing the
mutual influence among the cells’ resource allocations [14]–
[36]. The model, named load-coupling, refers to the timefrequency resource consumption in each cell is referred to
as the cell load. However, the model does not allow SIC. In
our recent work [1], we addressed resource optimization in
multi-cell NOMA. However, the system model is constrained
by fixed power allocation. How to model joint optimization
of power allocation and user pairing and how to solve the
resulting problems to optimality have remained open so far.
II. M AIN R ESULTS
We demonstrate how NOMA can be modeled in multi-cell
scenarios by significantly extending the approaches in [14]–
[36], with joint optimization of power allocation and user
pairing. One fundamental result under such type of models
in OMA is the existence of the equilibrium for resource allocation. However, this modeling approach in NOMA leads to
non-closed form formulation for cell load coupling, unlike the
case of OMA. The fact poses significant challenges in analysis
and problem solving. As one of our main results, we prove that
such an equilibrium for resource allocation in NOMA exists
as well and propose an efficient algorithm for obtaining the
equilibrium. Furthermore, we prove that the equilibrium is the
global optimum for resource optimization in multi-cell NOMA
and thus a wide class of resource optimization problems can be
optimally solved by our algorithm. Because of our analytical
results, previous works about OMA with load coupling is
a special case of ours, namely, the algorithmic notions and
mathematical tools being used in those works of classic multicell power control or OMA load coupling thus directly apply
to the analysis multi-cell NOMA, suggesting future works on
this topic.
To the best of our knowledge, this is the first work investigating how to optimally utilize power and time-frequency
resources jointly in multi-cell NOMA. As a key strength of
our modeling approach, it enables to formulate and optimize
an entire class of resource optimization problems. Namely,

2

III. C ELL L OAD M ODELING
Denote by C = {1, 2, . . . , n} and J = {1, 2, . . . , m} the sets
of cells and user equipments (UEs), respectively. Denote by
Ji the set of UEs served by cell i (i ∈ C). When using j to
refer to one UE in J, i by default indicates j’s serving cell,
unless stated otherwise. Downlink is considered in our model.
A. Resource Sharing in NOMA
The resource in time-frequency domain is divided into resource blocks (RBs). In OMA, one RB can be accessed by only
one UE. In NOMA, multiple UEs can be clustered together to
access the same RB by SIC. Increasing the number of UEs in
SIC, however, leads to fast growing decoding complexity [2],
[3]. In previous works, it has been demonstrated that most
of the possible performance improvement by SIC is reached
by pairing as few as two UEs [2]–[7]. Pairing two UEs for
resource sharing is illustrated in Figure 1. UEs within one
pair share the same RB and the RBs allocated to different
pairs do not overlap. We use u as a generic notation for a user
pair (referred to as “pair” for simplicity). For cell i (i ∈ C),
denote by Ui the set of candidate pairs. Suppose
 there are in
total mi UEs in cell i. Then |Ui | is up to m2i . Denote
by Vj
S
(j ∈ J) the set of pairs
containing
UE
j.
Let
U
=
U
i∈C i (or
S
equivalently U = j∈J Vj ) be the set of candidate pairs of all
cells. Let s = |U|. If there is a need to differentiate between
pairs, we put indices on u, i.e., U = {u1 , u2 , u3 , . . . us }. Finally,
in our model, for UEs we allow for both OMA and NOMA

RB
b

9
b

b

4
1

10

b

b

8
5

b

7
b
b

b

b

3

Frequency

as long as the cost function is monotonically increasing in
the cells’ time-frequency resource consumption, our proposed
framework in multi-cell NOMA applies. Specifically, for solving this class of problems optimally, we derive a polynomialtime algorithm S-C ELL that gives the optimal power allocation
and user pairing, for any given input of inter-cell interference.
To address the dynamic coupling of inter-cell interference, we
derive a unified algorithmic framework M-C ELL that solves
the multi-cell resource optimization problems optimally. The
algorithm S-C ELL serves as a sub-routine and is iteratively
called by M-C ELL. We demonstrate theoretically the linear
convergence of this process.
The fundamental differences between our investigated problems and single-cell NOMA are summarized as follows. For
multiple cells, the resource allocation in one cell affects the
interference that the cell generates to other cells. The amounts
of required resource to meet the demand for all cells are
coupled together, rather than being independent to each other.
Optimizing resource allocation within one cell leads to a chain
reaction among all other cells. Individual optimization for the
cells results in sub-optimality and very inaccurate performance
analysis. For multi-cell NOMA, not only the time-frequency
resource allocation but also the power splits and user pairings
in all cells are coupled together for the same reason. Therefore,
joint optimization in NOMA leads to a rather complex problem
for analysis.
By numerical experiments, optimizing resource utilization
by our algorithm enlightens how much we can gain from
NOMA in terms of resource efficiency and load balancing.

2

6
Time

User Pairing

Resource Allocation

Figure 1. This figure illustrates user pairing and time-frequency resource
sharing. There are 10 UEs in one cell. Eight form four user pairs {1, 5},
{2, 7}, {3, 6}, {8, 9}, and the other two UEs 4 and 10 are unpaired. The
UEs within one pair share the same time-frequency resource as indicated by
the colors.

with SIC. For each UE, that which mode is used (or both can
be used) is determined by optimization. In the following, we
refer to these two modes as orthogonal RB allocation and nonorthogonal RB allocation, respectively. In general NOMA, we
include both modes.
B. NOMA Downlink
We first consider orthogonal RB allocation in NOMA. Let
pi be the transmission power per RB in cell i (i ∈ C). Denote
by gij the channel coefficient from cell i to UE j. The signalto-interference-and-noise ratio (SINR) is:
pi gij
.
(1)
γj = P
2
k∈C\{i} Ikj + σ
The term Ikj denotes the inter-cell interference from cell k to
UE j, and is possibly zero. This generic notation is used for the
sake of presentation. Later, we use the load-coupling model,
where the cell load that reflects the usage of RBs governs the
amount of interference. The term σ2 is the noise power.
We then consider non-orthogonal RB allocation in NOMA.
In [37] (Chapter 6.2.2, pp. 238) it is shown that, with superposition coding, one UE of pair u (u ∈ U) can decode the other
by SIC. When there is need to consider the decoding order
in u, to be intuitive, we use ⊕ to denote the UE that applies
interference cancellation, followed by decoding its own signal.
And denotes the UE that only decodes its own signal. Note
that both ⊕ and are generic notations and refer to the two
different users in any pair u (u ∈ U) in consideration. For any
pair u, pi is divided to q⊕u and q u (q⊕u + q u = pi ), with
q⊕u and q u being allocated to ⊕ and , respectively. (The
generic notation qju (j ∈ u) denotes the power allocated to UE
j.) We remark that ⊕ decodes ’s signal first and hence ’s
signal does not compose the interference for ⊕. The SINR of
⊕ is computed by (2).
q⊕u gi⊕
.
(2)
γ⊕u = P
2
k∈C\{i} Ik⊕ + σ
The UE

is subject to intra-cell interference from ⊕, i.e.,
q u gi
X
γ u=
.
(3)
q⊕u gi +
Ik +σ2
| {z }
k∈C\{i}
intra-cell
| {z }
inter-cell

3

Denote by q the power allocation of all candidate pairs:


q
q⊕u2 · · · q⊕us
q = ⊕u1
.
q u1 q u2 · · · q us

We use qu to represent the column of pair u (u ∈ U) in q,
named power split for u. We remark that it is not necessary to
use all the pairs in U for resource sharing. Whether or not a
pair would be put in use and allocated with RBs is determined
by optimization, discussed later in Section III-E. In addition,
we remark that the decoding order is not constrained by the
power split [37], even though by our numerical results, more
power is always allocated to in optimal solutions. The issue
of the influence of inter-cell interference on the decoding order
is addressed later in Section III-D.
C. Inter-cell Interference Modeling
The basic idea is to use the cells’ RB consumption levels to
characterize respectively the cell’s likelihood of interfering to
the others. The approach is specified as follows. Denote by ρk
the proportion of RBs allocated for serving UEs in cell k. The
intuition behind the model is partially explained by the two
extreme cases ρk = 1 and ρk = 0. If cell k is fully loaded,
meaning that all RBs are allocated, then ρk = 1. In the other
extreme case, cell k is idle and accordingly ρk = 0. Consider
any UE j served by cell i. The interference j receives from
cell k is Ikj = pk gkj or Ikj = 0 in the two cases, respectively.
In general, ρk serves as a scaling parameter for interference,
see (4). By the interference modeling approach, the cell load
directly translates to the scaling effect of interference and
therefore the same notation is used for both.
Ikj = pk gkj ρk .

(4)

Intuitively, ρk reflects the likelihood that a UE outside cell
k receives interference from k. Note that ρk in fact is the
amount of time-frequency resource consumption of cell k and
hence is referred to as the load of cell k.
We remark that this type of interference modeling approach
is a suitable approximation for network-level performance
analysis, which enables study of inter-cell interference in
large-scale multi-cell networks without having to modeling
micro-level interference. Detailed system-level simulations
(e.g. [24] and [32]) have shown that this type of modeling has
sufficient accuracy for cell-level interference characterization.
This approach has been widely used and is getting increasingly
popular [14]–[36], which however, to our best knowledge, are
limited to OMA. We provide analytical results in order to
extend the modeling approach to NOMA.

for h, followed by decoding its own signal, and, user h does
not apply SIC. That is, the decoding order ⊕ = j and = h
always hold for pair independent of interference.
Proof. Denote by γhj and γhh respectively in (5) and (6) the
SINRs at users j and h for the downlink signal of h.
q g
P hu ij
.
(5)
γhj =
qju gij + k∈C\{i} pk gkj ρk + σ2
γhh =

qju gih +

Lemma 1. For u = {j, h} (gij > gih ) in cell i, if gij /gih >
gkj /gkh (k ∈ C\{i}), then SIC at j decodes first the signal

qhu gih
.
2
k∈C\{i} pk gkh ρk + σ

(6)

With superposition coding, j cancels the interference from h if
j can decode any data that h can decode [37], i.e. γhj > γhh ,
which reads:
X
qju gij gih + gij
pk gkh ρk + gij σ2
k∈C\{i}

> qju gij gih + gih

X

pk gkj ρk + gih σ2 .

k∈C\{i}

Further, γhj > γhh if and only if:
X pk ρk
(gih gkj − gij gkh ) 6 (gij − gih ).
σ2

(7)

k∈C\{i}

Recall that gij > gih , and therefore the right-hand side of (7)
is non-negative. Because of the condition gij /gih > gkj /gkh
for all k ∈ C\{i} in the statement of the lemma, the left-hand
side is non-positive. Hence Lemma 1.
The result of Lemma 1 is coherent with the previous
observations that two UEs with large difference in channel
conditions are preferred to be paired [4], [7]. If gij  gih ,
then most likely the condition in Lemma 1 holds, as the large
scale path-loss from other cells, tends not to differ as much
as from the serving cell i in this case. Besides, the large scale
path-loss is a practically reasonable factor for ranking the
decoding order [38], [39]. In Section VII, numerical results
further show that considering the UE pairs as defined by
Lemma 1 virtually does not lead to any loss of performance.
Lemma 1 is used to filter the candidate pairs set U (i.e. to drop
some candidate pairs from U) so as to reduce computational
complexity. From now on, we let Ui be composed of pairs
satisfying Lemma 1.
E. RB Allocation
If UE j (j ∈ J) is using orthogonal RB allocation, then the
achievable capacity1 of j is (8), with γj being (1).
cj = log(1 + γj ).

D. Determining Decoding Order
Inter-cell interference affects the decoding order in NOMA,
and thus how to model the load-coupling in NOMA is significantly more challenging than OMA. Lemma 1 below resolves
this issue by identifying pairs for which the decoding order
can be determined independently of interference. As another
benefit, it significantly reduces the set of candidate pairs.

P

(8)

For non-orthogonal RB allocation, the achievable capacity
for j and u (j ∈ u) is computed by cju = MB log (1 + γju )
with γju being (2) or (3). Therefore,
cju =

log (1 + γju )
0

j∈u
.
j∈
/u

(9)

1 For the sake of presentation, we use the natural logarithm throughout the
paper. We remark that all conclusions hold for the logarithm to base 2.

4

For UE j (j ∈ J), we use xj to denote the proportion of RBs
with orthogonal RB allocation to j. For any pair u (u ∈ U),
denote by xu the non-orthogonal RB allocation for the two UEs
in pair u. We use the vector x to represent the RB allocation
for all the UEs, i.e.,
x = [ x1 , x2 , . . . , xm ,
|
{z
}

xu , xu , . . . , xus ].
| 1 2{z
}

one can observe that increasing xu (or xj ) for some pair u
(or some UE j) may enhance the throughput of the UEs of u
(or UE j). However, since xu (or xj ) appears in the inter-cell
interference term (see (4) and (10)), the increase of xu (or xj )
results in less available resources for other UEs and leads to
more interference. The user pairing selection is not given a
priori but is determined by optimization.

Orthogonal RB allocation Non-orthogonal RB allocation

For any UE j, xj = 0 means that UE j does not use orthogonal
RB allocation. Similarly, for any pair u, xu = 0 means that
pair u is not put in use. For any UE j, if xu = 0 for all
u ∈ Vj , then it means that UE j only uses orthogonal RB
allocation. Resources used by different pairs are orthogonal
such that there is no interference among pairs. Denote by ρ̄
the cell load limit. By constraining that the sum of them which
equals to the load of cell i does not exceed ρ̄, the amounts
represented by xj (j ∈ Ji ) and xu (u ∈ Ui ) do not overlap.
Orthogonal RB allocation is considered among the pairs in
one cell, meaning that the pairs do not have interference with
each other.
X
X
xu 6 ρ̄.
(10)
ρi =
xj +
↑
↑
u∈Ui
j∈Ji
| {z }
| {z }
Cell
Load
load

Orthogonal
RB proportion

Non-orthogonal
RB proportion

ρ = [ρ1 , ρ2 , . . . , ρn ].
Similarly, we use vector ρ̄ to denote the load limits of all cells.
The term cj xj computes the bits delivered to UE j with
orthogonal RB allocation, because cj is the achievable capacity
of UE j on all RBs and xj is the proportion of RBs with
orthogonal RB allocation. Similarly, the term cju xu is the bits
delivered to UE j by non-orthogonal RB allocation for pair u.
Denote by dj the bits demand of UE j. The quality-of-service
(QoS) requirement is:
X
cj xj +
cju xu > dj .
(11)
|{z}
u∈V
↑
j
Bits delivered
| {z }
by orthogonal
Bits
Bits delivered by
non-orthogonal
RB allocation

The models proposed for OMA in [14]–[20] are inherently
a special case of our NOMA model. The former is obtained by
setting U = φ. Then, the terms for non-orthogonal RB allocation disappear in (10) and (11) and x is therefore eliminated
in (8)–(11). Also, there is no power split in OMA. Hence (8)–
(11) form a non-linear system only in terms of ρ. This system
falls into the analytical framework of standard interference
function (SIF) [40], which enables the computation of the
optimal network load settings via fixed-point iterations [17].
However, for the general NOMA case, the resource allocation
is not at UE-level. One needs to split a UE’s demand between
orthogonal and non-orthogonal RB allocations, which results
in a new dimension of complexity.

limit

We use ρ to represent the vector of network load, i.e.,

RB allocation

F. Comparison to OMA Modeling

IV. P ROBLEM F ORMULATION
By successively plugging (1) and (4) into (8), we obtain a
function cj in load ρ, i.e., cj (ρ). Similarly, we obtain cju (q, ρ)
from (2), (3), (4), and (9). For pair u (u ∈ U), we use a binary
variable yu to indicate whether or not the pair u is selected.
Define y as
y = [yu1 , yu2 , . . . , yus ].
We minimize a generic cost function F(ρ) that is monotonically (but not necessarily strictly monotonically) increasing in
each element of ρ. M IN F is given below.
[M IN F]
s.t.

min F(ρ)

ρi 6 ρ̄, i ∈ C
X
X
ρi =
xj +
xu , i ∈ C

demand

We remark that dj is normalized by the RB spectral bandwidth
and the total number of RBs, for the sake of presentation. Note
that a user can use orthogonal RB allocation individually, or
non-orthogonal RB allocation with the other user in the pair, or
both, which is subject to optimization. The amount of allocated
RBs to a user in OMA or a pair adopting NOMA, is subject
to optimization, under the constraint that the overall allocated
resource does not exceed limit.
We remark that there is an implicit pair selection problem
in the above expressions. Note that |U| increases fast with |J|.
It is therefore impractical to simultaneously use all pairs in U.
To deal with this issue, each UE is allowed to use up to one
pair in U for optimization, as formulated later in Section IV,
though our system model is not limited by this. The problem
of pairing and resource allocation is challenging: First, UEs
of the same pair are coupled in resource allocation. Second,

(12a)

ρ,q,x,y

j∈Ji

cj (ρ)xj +
X

X

(12b)
(12c)

u∈Ui

cju (q, ρ)xu > dj , j ∈ J

(12d)

u∈Vj

qju = pi , u ∈ Ui , i ∈ C

(12e)

j∈u

xu 6 yu , u ∈ U
X
yu 6 1, j ∈ J

(12f)
(12g)

u∈Vj

ρ, q, x > 0

(12h)

yu ∈ {0, 1}, u ∈ U

(12i)

Constraints (12b) guarantee that the cell load complies
to the load limit ρ̄. Constraints (12c) state the relationship
between RB allocation and cell load. Constraints (12d) and
(12e) are for QoS and power, respectively. Constraints (12f)
guarantee that RB allocation occurs only for selected pairs.
By constraints (12g), each UE belongs up to one pair such

5

that the selected pairs are mutually exclusive. The variables
are cell load ρ, power allocation q, RB allocation x, and
pair selection y. The variable domains are imposed by (12h)
and (12i). Throughout this paper, we use 0 to represent zero
vector/matrix. For simplicity, the dimension(s) of 0 is not
explicitly stated.
V. O PTIMIZATION WITHIN A C ELL
In multi-cell NOMA, due to the interference among cells,
one cell’s pair selection may affect the other cells’ power
splits, and vice versa. Let us consider a simple case in this
section. Suppose we optimize the load of one cell i, and the
cell load levels of C\{i} are temporarily fixed. This optimization step is a module for solving M IN F later in Section VI.
We respectively use qi , xi , yi to denote the corresponding
variable elements for power allocation, RB allocation, and pair
selection. Vector ρ−i is composed of all elements but ρi of
ρ. We minimize ρi under fixed ρ−i , as formulated below.
min

ρi ,qi ,xi ,yi

ρi s.t. (12c)–(12i) of cell i, with fixed ρ−i . (13)

Since ρ−i is fixed, cj is a constant and cju is a function in qi
only.
The optimization is not straightforward even under fixed
inter-cell interference. The optimal power split for one pair
is up to how much time-frequency resource is allocated to
this pair. In other words, for one pair u, if the amount of
RBs allocated to u changes, the optimal power split for u
before this change loses its optimality. So the power split q
and the resource allocation x are coupled together. In addition,
the pair selection is a combinatorial problem. Therefore, the
power split q, the time-frequency resource allocation x, and
the user pair selection y, must be optimized jointly.
Lemma 2. All constraints of (12d) in (13) hold as equalities
at any optimum.
Proof. Denote the optimal objective value of (13) by ρi0 and
the optimal orthogonal RB allocation of j (j ∈ Ji ) by xj0 .
Suppose strict inequality holds for some j. If xj0 > 0, by fixing
all other variables except for xj in (13), one can verify that
the solution xj0 −  ( > 0) is still feasible to (13) as long
as  is sufficiently small. In addition, xj0 −  leads to a lower
objective value ρi0 − , which contradicts our assumption that
ρi0 is the optimal objective value. If xj0 = 0, then j’s demand
has to be satisfied by non-orthogonal RB allocation and the
same argument applies to variable xu (j ∈ u).
The first analytical result is that, the optimal power split is
independent of pair selection, in Theorem 1. Denote by Yu the
set of all possible pairing solutions of (13) that includes pair
u, i.e.,
X
Yu = {yi |yu = 1,
yu 6 1, j ∈ Ji }.
u∈Vj

^ i (y
^ i ∈ Yu ), the optimal
Definition 1. Given a pair selection y
^ u , is the column
power split for pair u (u ∈ Ui ), denoted by q
^ i , where q
^ i is obtained by optimally solving
for pair u in q
^ i.
(13) for y

Theorem 1. Consider u (u ∈ Ui ). The optimal power split
^ i (y
^ i ∈ Yu ) is also optimal for y
^ i0 (y
^ i0 ∈ Yu ).
for any y
^ u and q
^ u0 the optimal power splits for
Proof. Denote by q
0
^ i and y
^ i , respectively. Suppose q
^ u is not optimal for y
^ i0 .
y
0
^ u ) = cju (q
^ u ) (j ∈ u); 2)
There are two possibilities: 1) cju (q
^ u ) 6= cju (q
^ u0 ) for at least one j in u.
cju (q
^ u and q
^ u0 result in the same xj and xu for
For 1), q
satisfying (12d) and are equally good for (13), which conflicts
^ u is optimal for y
^ i0 . We then consider
our assumption. Thus q
0
^ u ) > cju (q
^ u ). By Lemma 2, q
^ u0 makes
2) and assume cju (q
0
0
^ i . Replacing q
^ u by q
^ u leads to
(12d) become equality under y
some slack in (12d) and hence the objective can be improved.
^ u0 is optimal for yi0 . The same proof
This contradicts that q
^ u ) < cju (q
^ u0 ). Hence the conclusion.
applies to cju (q
By Theorem 1, the optimal power split is decoupled from
pair selection. Next we analytically prove how to find the
optimal power split for any pair.
A. Finding Optimal Power Split
Under fixed yi (yi ∈ Yu , u ∈ Ui ), constraints (12g)
are removed. Constraints (12f), and (12i) of (13) for all u
with yu = 0 in yi are removed. Therefore, for each pair
u = {⊕, }, we can formulate a problem in (14). Solving this
problem yields the optimal power split. In (14), x⊕ and x
are the orthogonal RB allocation for ⊕ and , respectively.
The variable xu denotes the amount of non-orthogonal RB
allocation for u.
x⊕ + x + xu

(14a)

c⊕ (ρ−i )x⊕ + c⊕u (qu , ρ−i )xu > d⊕

(14b)

min

x⊕ ,x ,xu >0
qu >0

s.t.

c (ρ−i )x + c
q⊕u + q

u

u (qu , ρ−i )xu

>d

= pi

(14c)
(14d)

For deriving solution method for (14), define function wj
(j = ⊕ or j = ) of ρ−i as follows.

,
X
wj (ρ−i ) = 
pk gkj ρk + σ2  gij .
(15)
k∈C\{i}

For q⊕u , one can derive from (2) and (9):

q⊕u = (ec⊕u − 1)w⊕ (ρ−i ).

(16)

Combining (16) with (14d), q⊕u and q u can be eliminated,
giving (17) below. Formulation (17) is equivalent to (14).
Given c⊕u and c u , the corresponding q⊕u and q u can be
obtained from c⊕u and c u by (14d) and (16).
min

x⊕ ,x ,xu >0
c⊕u ,c u >0

x⊕ + x + xu

s.t. c⊕ (ρ−i )x⊕ + c⊕u xu > d⊕
c (ρ−i )x + c
Cvu (c⊕u , c

u xu > d

u , ρ−i )

60

(17a)
(17b)
(17c)
(17d)

6

In (17), the function Cvu is defined in (32) in the Appendix.
One can easily verify that Cvu is convex in c⊕u and c u (with
gi⊕ > gi ). The difficulty of (17) is on the two bi-linear
constraints (17b) and (17c). However, they become linear with
fixed xu . To ease the presentation, we define the function
below.
Zu (xu , ρ−i ) = xu +

min

x⊕ ,x >0
c⊕u ,c u >0

x⊕ + x

c⊖u

K
K= [cK
⊕u , c⊖u ]

s.t. (17b)–(17d). (18)

Solving (17) (and equivalently (14)) is to find the minimum
of Zu (xu , ρ−i ). The following theorem shows the uniqueness
of the minimum of Zu (xu , ρ−i ).

[0, 0]

Theorem 2. Zu (xu , ρ−i ) has unique minimum in xu .
Proof. Since the first term xu in Zu (xu , ρ−i ) is strictly monotonically increasing, to prove that Zu (xu , ρ−i ) has unique
minimum, we only need to prove that the remaining part
of Zu (xu , ρ−i ) is monotonically (but not necessarily strictly
monotonically) decreasing in xu . For this part, at the optimum (17b) and (17c) hold as equalities because of Lemma 2.
Hence, reformulating the problem by replacing the inequalities
in (17b) and (17c) with equalities does not lose optimality.
With equalities, the variables x⊕ and x can be represented
by c⊕u and c u :
(d − c u xu )
(d⊕ − c⊕u xu )
, x =
.
x⊕ =
c⊕ (ρ−i )
c (ρ−i )

(19)

Therefore x⊕ and x can be eliminated from the objective
function. The minimization is thus equivalent to maximizing
c⊕u /c⊕ (ρ−i ) + c u /c (ρ−i ). We formulate this maximization problem below.
c⊕u
c u
max
+
(20a)
c⊕u ,c u >0 c⊕ (ρ−i )
c (ρ−i )
s.t. c⊕u xu 6 d⊕
(20b)
c

u xu

6d

Cvu (c⊕u , c

(20c)
u , ρ−i )

60

(20d)

Constraints (20b) and (20c) originate from the nonnegativity requirement of x⊕ and x . Note that (20) is convex.
In addition, the feasible region shrinks with the increase of xu .
Then the optimum of (20) monotonically decreases with xu .
Hence the theorem.
Because of Theorem 2, bi-section search of xu reaches the
minimum of Zu (xu , ρ−i ). Note that for any xu , computing
Zu (xu , ρ−i ) needs to solve (20). In the following, we prove
how this can be done much more efficiently than employing
standard convex optimization.
We remark that (20) is a two-dimensional optimization
problem with respect to c⊕u and c u . Constraints (20b)
and (20c) are defined by two hyperplanes c⊕u = d⊕ /xu and
c u = d /xu , respectively. Due to the convexity of the function Cvu in c⊕u and c u , the curve Cvu (c⊕u , c u , ρ−i ) = 0
in (20d) along with c⊕u = 0 and c u = 0 forms a convex
set of [c⊕u , c u ]. The optimum of (20) depends on whether
the two hyperplanes intersect with the curve and how they
intersect. This leads to three possible cases to be considered,

h

1
1
c⊕ (ρ−i ) , c⊖ (ρ−i )

i

c⊕u

Figure 2. This figure shows the feasible region of (20) of c⊕u and c u
for different xu . The shadowed area below the curve is (20d). The vertical
and horizontal dashed lines are the hyperplanes defined by (20b) and (20c),
respectively. The hyperplanes are shown by red, blue, and green dashed lines
for Case 1, Case 2, and Case 3, respectively. The point K is optimal for Case
1. The blue and green circles represent the optimal solutions for Case 2 and
Case 3, respectively. The black dot is the point on the curve where the two
hyperplanes intersect with each other (see footnote 3 for the existence of this
point.).

named Case 1, Case 2, and Case 3, respectively. In Case 1,
constraints (20b) and (20c) are redundant, and the optimum is
determined by the coefficients 1/c⊕ (ρ−i ) and 1/c (ρ−i ) in
the objective function and the curve Cvu (c⊕u , c u , ρ−i ) = 0.
In Case 2, the optimum is defined by one of the hyperplanes
and the curve. In Case 3, constraint (20d) is redundant, and
the optimum is determined by the two hyperplanes. With xu
increasing from 0 to ∞, Case 1, Case 2, and Case 3 happen
sequentially, and all happen eventually. The three cases are
illustrated in Figure 2 with colors. Below we respectively show
how to compute the optimum for each case.
In the following, we first compute the optimum in Case 1,
K
represented by K = [cK
⊕u , c u ], which is the intersection of the
vector [1/c⊕ (ρ−i ), 1/c (ρ−i )] and the curve. The point also
leads to a closed-form solution for the optima of all the three
cases. Mathematically, point K is solved by applying bi-section
search to (21) below.
Cvu (c⊕u , c u , ρ−i ) = 0
c⊕u c⊕ (ρ−i ) = c u c (ρ−i ).

(21)

For solving (21), we can first eliminate c⊕u or c u in the first
equation. This can be done by representing one of c⊕u and
c u with the other by the second equation. Since there is only
one variable in the first equation, one can use bi-section search
to find its solution2 . Then, we compute the value of xu when
at least one hyperplane goes through K, denoted by xK
u in (22).
K
K
xK
u = min{d⊕ /c⊕u , d /c u }.

(22)

The three cases, indicated by colors in Figure 2, are as follows.
2 The solution is guaranteed to be unique and hence bi-section search
applies. This is because, by representing one of c⊕u and c u by the other by
function Cvu , one variable is monotonically decreasing in the other, resulting
in a unique zero point.

7

Case 1 (xu 6 minj∈u dj /cK
ju ): Point K is the optimum of
(20), because (20b) and (20c) are redundant. This happens
when xu is sufficiently small (or 0), as shown in Figure 2.
Case 2 (xu > minj∈u dj /cK
ju and Cvu (d⊕ /xu , d /xu ) > 0):
There exists one point on the curve where both two hyperplanes intersect3 . We represent this point by the black dot
on the curve in Figure 2. In Case 2, one hyperplane intersects
with the curve at some point between K and the black dot, and
intersects with the other hyperplane on some point above the
curve, see Figure 2. The intersection point of the curve and the
hyperplane is the optimum of (20). Without loss of generality,
we assume K violates (20c), meaning that the hyperplane of
(20c) goes through the optimum, as shown by Figure 2. By
plugging the equation of the hyperplane into that of the curve,
the optimal c⊕u is a function of xu . Similarly, if K violates
(20c) instead, then c⊕u = d⊕ /xu and the optimal c u is
a function of xu . To know which hyperplane goes through
K
the optimum, one only needs to check which of cK
ju xju > dj
(j = ⊕ or j = ) holds. Note that exactly one of the two holds
in Case 2. The optimal c⊕u and c u are computed respectively
by (33) and (34) defined in the Appendix.
Case 3 (xu > minj∈u dj /cK
ju and Cvu (d⊕ /xu , d /xu ) 6 0):
Constraint (20d) is redundant, as shown in Figure 2. The
optimum is the intersection point of the two hyperplanes,
computed by cju = dj /xu (j = ⊕ or j = ).
In summary, the optimal solution of (20) is computed by
(23) below (j = ⊕ or j = ) in closed form, with Hju
being (33) or (34) in the Appendix.
 K
Case 1
 cju
Cju (xu , ρ−i ) =
(23)
H (x , ρ ) Case 2
 ju u −i
dj /xu
Case 3
The function Zu (xu , ρ−i ) computes the amount of resource
used for both orthogonal and non-orthogonal RB allocations
for the UEs in u. It is optimal to serve the two UEs only
by orthogonal RB allocation, if the minimum of Zu (xu , ρ−i )
occurs at xu = 0. In all other cases, minxu Zu (xu , ρ−i ) yields
the optimal power split for non-orthogonal RB allocations. The
algorithm optimally solving (14), named S PLIT, is as follows.
S PLIT(u, ρ−i )
1 x∗u = arg minxu Zu (xu , ρ−i ) // Bi-section search
2 Compute hc∗⊕u , c∗ u i by (23) // With x∗u , ρ−i
3 Compute hx∗⊕ , x∗ i by (19) // With c∗⊕u , c∗ u , ρ−i
4 Convert hc∗⊕u , c∗ u i to hq∗⊕u , q∗ u i // By (14d), (16)
5 return hq∗⊕u , q∗ u , x∗⊕ , x∗ , x∗u i

B. Optimal Pairing
Denote by Yi the set of all candidate pair selections:
Yi = (∪u∈Ui Yu ) ∪ {0}.
3 The existence of this point is guaranteed: With the increase of x , both
u
hyperplanes will eventually intersect with the curve with two intersection
points. By increasing xu , the distance between the two intersections keeps
being smaller. The two intersections will eventually overlap.

0.18
1
0.02

∆

0.14

0.14
0.16 3

2
0.07

0.06

0.40
4

5
0.12

Figure 3. The figure shows an example of one cell i with five UEs,
i.e., Ji = {1, 2, 3, 4, 5}. Assume the candidate pair set is Ui =
{{1, 2}, {2, 4}, {4, 5}, {3, 5}}. The blue edges are type-1. The red edges
are type-2. A matching is a set of edges without common vertices (also called
independent edge set) and is a pair selection solution. The maximum matching,
as highlighted in the figure, is {{1, 2}, {3, 5}, {4, ∆}}.

By obtaining minxu Zu (xu , ρ−i ) for all u ∈ Ui as shown
earlier in Section V-A, enumerating all yi in Yi gives the
optimal solution to (13). This exhaustive search however does
not scale, as |Yi | is exponential in the number of UEs. By
the following derivation, we are able to obtain the optimum
of (13) in polynomial time.
Theorem 3. The optimum of (13) is computed by finding the
maximum weighted matching in an undirected graph.
Proof. To prove the conclusion, an undirected weighted graph
Gi is constructed and explained below.
Gi =

hJi , Ui , wi
hJi ∪ {∆}, Ui ∪ {{j, ∆}|j ∈ Ji }, wi

|Ji | even
|Ji | odd.

(24)

In (24), the graph is represented by a 3-tuple, with the first
element being the vertex set, the second element being the
edge set, and the third element being the weight vector.
Parameter ∆ is an auxiliary vertex for odd |Ji |. Without loss
of generality, below we focus on odd |Ji |. (All conclusions
naturally hold for |Ji | being even.) By the definition in (24),
each UE is corresponding to a vertex. For each pair u in Ui ,
there is one edge connecting the two UEs in u, associated
with weight wu . We name these as type-1 edges. Besides, for
each UE j in Ji , there is one extra edge connecting j and the
auxiliary vertex ∆, associated with weight wj . We name these
as type-2 edges. An illustration is given in Figure 3.
The weight w is defined as follows, where T is a positive
value keeping all weights being positive.
Type-1 edge
Type-2 edge

wu = T − minxu Zu (xu , ρ−i ) (u ∈ Ui )
wj = T − dj /cj (ρ−i ) (j ∈ Ji )

First, we remark that any yi is feasible to (13) if and only
if all the pairs u with yu = 1 (u ∈ Ui ) form a matching (or
an
P empty edge set) in Gi . Otherwise, there exists j such that
u∈Vj yu > 2, and (12g) would be violated. Then, by the
definition of weights, minimizing the load ρi becomes finding
a maximum weighted matching.
The algorithm S-C ELL solving (13) exactly is as follows.
Lines 1–8 compute the edge weights of the graph to be

8

S-C ELL(ρ−i )
1 for u ∈ Ui
2
hq⊕u , q u , x⊕ , x , xu i = S PLIT(u, ρ−i )
3
wu = T − (x⊕ + x + xu ) // wu > 0
4 end for
5 if |Ji | is odd
6
for j ∈ Ji
7
xj = dj /cj (ρ−i )
8
wj = T − xj // wj > 0
9
end for
10 end if
11 Construct Gi by (24)
12 U∗i = M AXIMUM -W EIGHTED -M ATCHING(Gi )
13 x∗i = 0, q∗i = 0, y∗i = 0
14 for u ∈ U∗i ∩ Ui
15
x∗u = xu , y∗u = 1
16
for j ∈ u
17
q∗ju = qju , x∗j = xj
18
end for
19 end for
20 if |Ji | is odd
21
Find the {j, ∆} in U∗i and let x∗j = xj
22 end ifP
P
23 ρi = j∈Ji x∗j + u∈Ui x∗u
24 return hρi , q∗i , x∗i , y∗i i

To prove (13) is feasible, we prove the remaining problem is
always feasible. Note that, with yi and qu being fixed, (13)
becomes a linear programming (LP) problem, which is stated
below (the equalities are by Lemma 2).
X
X
X
min
xj +
xu , s.t. cj xj +
cju xu = dj . (25)
x>0

u∈Ui

j∈Ji

u∈Vj

By Farkas’ lemma, a group of linear constraints in standard
form, i.e. Ax = b (b > 0), is feasible with x > 0 if and only
if there does not exist v such that v> A > 0> and v> b < 0.
Obviously, there is no v with d > 0 satisfying (26).
X
vj > 0 (j ∈ Ji ) and
v j dj < 0
(26)
j∈Ji

Hence (25) is feasible, and the conclusion holds.
Let λi = |Yi |. For each yi in Yi , we use an integer in
[1, λi ] to uniquely index yi . We refer to all the pair selection
solutions in Yi as pairing 1, pairing 2, . . . , pairing λi . Denote
by fik (ρ−i ) the optimum of (13) under pairing k (1 6 k 6
λi ), i.e.,
fik (ρ−i ) =

min ρi s.t. (12c)–(12f) and (12h) of cell i.

ρi ,qi ,xi

Let fi (ρ−i ) be the optimum5 of (13). Then we have:
fi (ρ−i ) =

min

k=1,2,...,λi

fik (ρ−i ).

Network-wisely, we have:
f(ρ) = [f1 (ρ−1 ), f2 (ρ−2 ), . . . , fn (ρ−n )].

constructed. Then we construct the graph in Line 11 and
compute the maximum matching4 U∗i in Line 12, which by
Theorem 3 is the optimal pair selection in cell i. Lines 13–22
assign the obtained solutions to hq∗⊕u , q∗ u , x∗⊕ , x∗ , x∗u i for the
pairs in U∗i . The other pairs are not selected and hence their
values in x∗i , q∗i , and y∗i are zeros.
We remark that in the matching process, if the number
of nodes is odd, for the unpaired UE, it is allocated with
orthogonal RBs. For two UEs that are paired in the solution of
matching, they are not necessarily in non-orthogonal allocation
but is up to optimization.
VI. M ULTI - CELL L OAD O PTIMIZATION
This section proposes the algorithmic framework M-C ELL
for deriving the optimum of M IN F, by analyzing sufficientand-necessary conditions of optimality and feasibility.
A. Revisiting Single-cell Load Minimization
Recall that for single cell optimization, the optimum of (13)
of cell i (i ∈ C) is a function of the load of other cells ρ−i .
By Lemma 3 below, this function is well-defined for any nonnegative ρ−i .

best known algorithm
[41] runs on Gi in O((|Ui | + |Ji |)
p
(odd |Ji |) or O(|Ui | |Ji |) (even |Ji |).
4 The

Theorem 4. f(ρ) is an SIF, i.e. the following properties hold:
1) (Scalability) αf(ρ) > f(αρ), ρ > 0, α > 1.
2) (Monotonicity) f(ρ) > f(ρ 0 ), ρ > ρ 0 , ρ, ρ 0 > 0.
Proof. We first prove monotonicity and scalability for
fik (ρ−i ) (i ∈ C, k = 1, 2, . . . , λi ). For monotonicity, we
0
0
prove that fik (ρ−i
) 6 fik (ρ−i ) for ρ−i
6 ρ−i as follows.
0
Given any non-negative ρ−i , we replace ρ−i with ρ−i
. Note
0
that cju (ρ−i ) > cju (ρ−i ). Thus the replacement makes the
0
solution space of (12d) larger, and the optimum with ρ−i
is no
0
larger than that with ρ−i . Therefore fik (ρ−i ) 6 fik (ρ−i ). For
scalability, we prove that fik (αρ−i ) 6 αfik (ρ−i ) for α > 1
and non-negative ρ−i as follows. Denote the optimal solution
of fik (ρ−i ) by hρi00 , qi00 , xi00 i. We have fik (ρ−i ) = ρi00 . Due
to that 1/cju (qi , ρ−i ) and 1/cj (ρ−i ) are strictly concave in
ρ−i , the two inequalities
α
1
α
1
<
,
<
(29)
cj (αρ−i )
cj (ρ−i ) cju (qi , αρ−i )
cju (qi , ρ−i )
hold for α > 1. Consider the following minimization problem 30, with yi being fixed to pairing k.
min

ρi ,qi ,xi >0

p
|Ji |)

(28)

The following theorem reveals a key property of f(ρ).

Lemma 3. The problem in (13) is always feasible.
Proof. We select some yi in Yi and fix it in (13) (Yi 6= φ
by definition). For each pair u, we fix qu to [pi /2, pi /2]> .

(27)

s.t.

ρi

(30a)

(12c), (12e) and (12f) of cell i
(30b)
X
cj (ρ−i )xj +
cju (q, ρ−i )xu > αdj , j ∈ Ji (30c)
u∈Vj

5 Therefore,

fi (ρ−i ) equals the ρi obtained from S-C ELL(ρ−i ).

9

Note that hαρi00 , qi00 , αxi00 i is feasible to (30), with the objective
value being αfik (ρ−i ). Hence the optimum of (30) is no more
than αfik (ρ−i ). For fik (αρ−i ), note that the corresponding
optimization problem only differs with (30) in (30c). Instead
of (30c), in fik (αρ−i ) we have:
X
cj (αρ−i )xj +
cju (q, αρ−i ) > dj , j ∈ Ji
(31)
u∈Vj

By Lemma 2, (30c) is equality at the optimum. Then by (29),
for any solution of (30), using it for the optimization problem
associated with fik (αρ−i ) makes (31) an inequality. Therefore the problem for fik (αρ−i ) has a lower optimum than
(30). Further, the optimum is lower than αfik (ρ−i ). Hence
fik (αρ−i ) < αfik (ρ−i ).
We then allow k to be variable and consider fi (ρ−i ) (i ∈ C).
0
For ρ−i
6 ρ−i we have
0
0
fi (ρ−i
) = min fik (ρ−i
) 6 min fik (ρ−i ) = fi (ρ−i )
k

k

and for α > 1 we have
fi (αρ−i ) = min fik (αρ−i ) < α min fik (ρ−i ) = αfi (ρ−i )
k

k

Hence the conclusion.
Given ρ, denote by fk (k > 1) the function composition of
f(fk−1 (ρ)) (with f0 (ρ) = ρ). Lemma 4 holds by [40].
Lemma 4. If limk→∞ fk (ρ) exists, then it exists uniquely for
any ρ > 0.
B. Optimality and Feasibility
Based on Theorem 4, we derive sufficient-and-necessary
conditions for M IN F in terms of its feasibility and optimality.
For any load ρ, we say that a load ρ is achievable if and only
if there exist q, x, and y such that the solution hρ, q, x, yi is
feasible to M IN F.
Lemma 5. For any ρ > 0, if there exists i ∈ C such that
ρi < fi (ρ−i ), then ρ is not achievable in M IN F.
Proof. Let ρi0 = fi (ρ−i ). By the definition of fi , ρi0 is the
minimum value satisfying (12c)–(12i) under ρ−i . Therefore
any ρi with ρi < ρi0 is not achievable with constraints (12c)–
(12i). Hence the conclusion.
Theorem 5. In M IN F, ρ (ρ 6 ρ̄) is achievable if and only if
f(ρ) is achievable and ρ > f(ρ).
Proof. By the inverse proposition of Lemma 5, an achievable
ρ always satisfies ρ > f(ρ). The necessity is proved as
follows. Suppose ρ is achievable for M IN F. Consider using
f(ρ) as another solution (together with the hq, x, yi obtained
when computing f(ρ)). Then f(ρ) satisfies (12b). Also, f(ρ)
together with its hq, x, yi fulfills (12c)–(12i) by the definition
of f(ρ). Thus, f(ρ) is achievable.
For the sufficiency, note that the achievability of f(ρ)
implies that ρ along with hq, x, yi obtained by solving f(ρ)
satisfies (12c)–(12i). Combined with the precondition ρi 6 ρ̄
(i ∈ C), the load ρ is feasible to (12b)–(12i) (and thus
achievable in M IN F). Hence the conclusion.

Theorem 5 provides an effective method for improving any
sub-optimal solution to M IN F. For any achievable ρ, evaluating f(ρ) always yields a better solution6 . This conclusion
is based on Theorem 5: Suppose ρ (ρ > 0) is the current
cell load, and let ρ 0 be the function value evaluated at ρ, i.e.
ρ 0 = f(ρ). By Theorem 5, we always have ρ 0 6 ρ.
Recall that F(ρ) is the objective function of the problem
M IN F. Theorem 6 below states that, the fixed point of f(ρ)
(along with hq, x, yi obtained when computing f(ρ)) is optimal to M IN F.
Theorem 6. Load ρ∗ is the optimum of M IN F if (and only if
when F(ρ) is strictly monotonic) ρ∗ = f(ρ∗ ) 6 ρ̄.
Proof. (Necessity) If ρ∗ is optimal (and thus feasible), then
obviously we have ρ∗ 6 ρ̄. By Theorem 5, f(ρ∗ ) is also
feasible and f(ρ∗ ) 6 ρ∗ . By successively applying Theorem 5,
fk (ρ∗ ) for any k > 1 is a feasible solution and fk (ρ∗ ) 6
fk−1 (ρ∗ ). Let ρ 0 = limk→∞ fk (ρ∗ ). Then ρ 0 6 ρ∗ holds by
the above derivation. In addition, note that ρ 0 is a feasible
solution as well. By that ρ∗ is optimal for M IN F, we have
ρ 0 = ρ∗ , otherwise ρ 0 would lead to a better objective value in
M IN F than ρ∗ . Hence ρ∗ = limk→∞ fk (ρ∗ ), i.e. ρ∗ = f(ρ∗ ).
(Sufficiency) By Theorem 5, for any feasible ρ,
limk→∞ fk (ρ) is feasible and limk→∞ fk (ρ) 6 ρ holds.
By Lemma 4, the limit remains for any ρ > 0, and thus
limk→∞ fk (ρ) = limk→∞ fk (ρ∗ ). Since ρ∗ = f(ρ∗ ), we
have ρ∗ = limk→∞ fk (ρ∗ ). Thus ρ∗ 6 ρ for any feasible
ρ, meaning that ρ∗ is optimal for M IN F.
Hence the conclusion.
C. The Algorithmic Framework
Starting from any non-negative ρ(0) , we compute
limk→∞ fk (ρ) iteratively. During each iteration, n problems
in (13) for i ∈ C are solved. The convergence is guaranteed by
Lemma 4. At the convergence, by Theorem 6, the optimum is
reached. Note that once ρ(k) is feasible for any k > 0, then by
Theorem 5, all ρ(k+1) , ρ(k+2) , . . . are feasible as well. One
can terminate prematurely to obtain a sub-optimal solution
with less computation. M-C ELL is outlined below.
M-C ELL(ρ(0) , )
1 k=0
2 repeat
3
k = k+1
4
for i ∈ C
(k)
(k) (k)
(k)
(k−1)
5
hρi , qi , xi , yi i = S-C ELL(ρ−i )
6
end for
7 until kρ(k) − ρ(k−1) k∞ 6 
(k)
8 if ρi > ρ̄ for some i (i ∈ C)
9
M IN F is infeasible
10 return hρ(k) , q(k) , x(k) , y(k) i
M-C ELL applies fixed point iterations using f(ρ). The
convergence of fixed point iterations on f(ρ) is linear [42]. The
6 Rigorously, Theorem 5 implies that the new solution is not worse. In fact
it is guaranteed to be strictly better (with strictly monotonic F(ρ)) unless the
old one is already optimal. A proof can be easily derived based on Theorem 6.

10

Total load

15

called “asynchronous fixed-point iterations” [40] can be used.
The asynchronous fixed-point iterations converge to the fixed
point that is the same as obtained by its synchronized version.
Intuitively, the fixed point is unique, regardless of how we
reach it.

OMA(Opt)
SP(Uni)
SP(FTPC)
NOMA(Opt)

10

VII. P ERFORMANCE E VALUATION

5

0

0.2

0.4

0.6

0.8

1

Normalized demand d

Max load

1

OMA(Opt)
SP(Uni)
SP(FTPC)
NOMA(Opt)

We use a cellular network of 19 cells. To eliminate edge
effects, wrap-around technique [43] is applied. Inside each
cell, 30 UEs are randomly and uniformly distributed. In each
cell, there are in total 30
= 435 possible choices for user
2
pairing in NOMA. User demands are set to be a uniform
value d. In the simulations, d is normalized by M × B in (8)
and (9), and belongs to (0, 1]. The network in OMA reaches
the resource limit at d = 1.0, i.e., any d > 1.0 leads to at
least one cell being overload in OMA. Other parameters are
given in Table I.
Table I
S IMULATION PARAMETERS .

0.5

0
0.2

0.4

0.6

0.8

1

Normalized demand d
Figure 4. This figure illustrates the total and maximum load in function of
normalized demand. At d = 1.0, the network reaches its resource limit
such that any larger demand cannot be satisfied by OMA(Opt) . SP(Uni) and
SP(FTPC) are two sub-optimal NOMA power allocation schemes, for which,
pair selections are optimally computed.

feasibility check is done by Lines 8 and 9. The infeasibility
of M IN F implies that at least one cell will be overloaded for
meeting user demands. If this happens, we know for sure
that the user demands cannot be satisfied. We remark that
all the conclusions derived in this section are independent
of the implementation of S-C ELL in Line 5. As long as the
sub-routine S-C ELL yields the optimal solution to (13), MC ELL achieves the optimum of M IN F7 . Besides, M-C ELL
possesses the optimality for M IN F with any objective function
that is monotonically (but not necessarily strictly monotonic)
increasing in each element of ρ. These two properties make
M-C ELL an algorithmic framework. To our knowledge, the
most efficient S-C ELL is what we derived in Section V.
For a cell i (i ∈ C), given the information of other cells’ load
ρ−i , solving fi (ρ) is based on local information, making MC ELL suitable to run in a distributed manner. A cell can maintain the information of a subset of cells (e.g., the surrounding
cells) having major significance in terms of interference, and
exchange the information with other cells periodically, which
can be implemented via the LTE X2 interface. The technique
7 With filtered U, the proposed M-C ELL is proved to converge to the global
optimum of M IN F. Without filtered U (or for any possible candidate pairs
set U), M-C ELL is still applicable to M IN F though there is no theoretical
guarantee of convergence or optimality, as the decoding order for each pair
may change in the iteration process.

Parameter
Cell radius
Carrier frequency
Total bandwidth
Cell load limit ρ̄
Path loss model
Shadowing (Log-normal)
Fading
Noise power spectral density
RB power pi (i ∈ C)
Convergence tolerance ()

Value
500 m
2 GHz
20 MHz
1.0
COST-231-HATA
6 dB standard deviation
Rayleigh flat fading
−173 dBm/Hz
800 mW
10−4

We consider two objectives for performance evaluation:
resource efficiency and load balancing.
P For resource efficiency,
the objective function is F(ρ) = i∈C ρi , i.e., to minimize
the total network time-frequency resource consumption (or
cells’ average resource consumption if divided by n). For load
balancing, we adopt min-max fairness and the objective function is F(ρ) = maxi∈C ρi . Section VII-A and Section VII-B
provide results for power allocation and user pairing, respectively. The optimal OMA, named OMA(Opt) , is obtained by
fixing y to 0 in M IN F and solving the remaining problem to
optimality8 . The proposed optimal NOMA solution is named
NOMA(Opt) in the remaining context.
A. Power Allocation
We use OMA(Opt) as baseline. As for NOMA, the pairing
candidate set U initially covers all pairs of UEs in each cell.
Then, those pairs not fulfilling Lemma 1 are dropped from
U. We then use M-C ELL to compute NOMA(Opt) . Besides the
optimal NOMA, we implement two other sub-optimal NOMA
power split schemes for comparison. One is named “SP(Uni) ”,
in which the power pi splits equally between q⊕u and q u for
any pair u = {⊕, } (u ∈ U). The other is “fractional transmit
power control” (FTPC), named SP(FTPC) , using a parameter
8 With y being fixed to 0 in M IN F, the variables q and x disappear. Then
P
(k)
we modify Line 5 of M-C ELL to be “ρi =
j∈Ji dj /cj (ρ−i )” and
Line 10 to be “return ρ(k) ”. The modified M-C ELL gives the optimal load
for OMA (see [17] for further details).

11

Best

Best

(a) PA(B-W)

Best

Best

(c) Optimal pairing in filtered U

(b) PA(B-SB)

(d) Pairing in non-filtered U

Figure 5. This figure illustrates pair selection in a typical cell of 30 UEs. The UEs are represented by the vertices on a circle. The UE marked “Best” at the
top position has the best channel condition. The UEs are arranged clock-wisely in the descending order of channel conditions. The edges are selected pairs.
Each subfigure represents one pairing method. Figure 5(a) and Figure 5(b) show PA(B-W) and PA(B-SB) , respectively. Figure 5(c) shows optimal pairing after
filtered U. Figure 5(d) shows the pairing solution obtained with non-filtered U.

Cell load

0.7
0.6

|U| =

30
2 × 19

= 8265

|U| = 5779 (U filtered by Lemma 1)

0.5
0.4
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

Cell index
Figure 6. This figure shows the optimized
 load levels of all 19 cells. The cells are numbered in an ascending order of loads. The blue bars show the computed
cell load under U composed of all 30
2 × 19 pairs. The red bars show the minimum cell load with pairs satisfying Lemma 1.

0.6
0.4

0.8

OMA(Opt)
PA(B-W) -SP(Uni)
PA(B-W) -SP(FTPC)
PA(B-W) -SP(Opt)
NOMA(Opt)

Average load

Average load

0.8

0.2

0.6
0.4

OMA(Opt)
PA(B-SB) -SP(Uni)
PA(B-SB) -SP(FTPC)
PA(B-SB) -SP(Opt)
NOMA(Opt)

0.2

0

0
0.2

0.4

0.6

0.8

1.0

Normalized demand d

0.2

0.4

0.6

0.8

1.0

Normalized demand d

Figure 7. This figure evaluates PA(B-W) , combined with three power split
schemes, SP(Uni) , SP(FTPC) , and SP(Opt) . In SP(Opt) , the power split is optimal
for each pair. OMA(Opt) and NOMA(Opt) are baselines.

Figure 8. This figure evaluates PA(B-SB) , combined with three power split
schemes, SP(Uni) , SP(FTPC) , and SP(Opt) . In SP(Opt) , the power split is optimal
for each pair. OMA(Opt) and NOMA(Opt) are baselines.

to control the fairness for power split. We set this parameter
to be 0.4 as recommended in [44]. Under both SP(Uni) and
SP(FTPC) , we use the method in Section V-B to compute the
optimal pair selection. Both two power split schemes are easily
accommodated by M-C ELL.
Figure 4 shows the total load and the maximum load in
function of normalized demand. As expected, the cell load
levels monotonically increase with user demand. At high user
demand, NOMA(Opt) dramatically improves the load performance. For d = 1.0, it achieves 31% better performance

than OMA(Opt) for both total load and maximum load. The
two sub-optimal solutions SP(Uni) and SP(FTPC) also result in
load improvement than OMA(Opt) . Compared to the two suboptimal solutions, the improvement achieved by NOMA(Opt)
over OMA(Opt) is doubled or more. On average, by using the
same amount of time-frequency resource, NOMA(Opt) delivers
33% more bits demand than OMA(Opt) . Besides, SP(FTPC)
achieves better performance than SP(Uni) , as the former takes
into account the channel conditions in power split. Generally,
in SP(FTPC) , UE with worse channel is allocated with more

12

B. User Pairing
We study the influence of user pairing by considering
two sub-optimal ones [4], named “PA(B-W) ” and “PA(B-SB) ”,
respectively. Suppose we sort the UEs in descending order
of their channel conditions. In PA(B-W) , the UE with the best
channel condition is paired with the UE with the worst, and
the UE with the second best is paired with one with the
second worst, and so on. In PA(B-SB) , the UE with the best
channel condition is paired with the one with the second
best, and so on. See Figures 5(a) and 5(b) for an illustration.
In addition, we examine to what extend pair filtering (by
Lemma 1) affects performance. For filtered U, optimal pair
selection is done by Section V-B. For non-filtered U, we
apply M-C ELL even though there is no theoretical guarantee
on optimality. Convergence, however, is observed for all the
instances we considered. Figures 5(c) and 5(d) illustrated the
resulted selection patterns.
In Figure 6, we show the load levels of all 19 cells with
d = 1.0, under both filtered and non-filtered
U. In this

specific scenario, |U| is reduced from 30
2 × 19 = 8265 to
5779 after being filtered by Lemma 1. We choose d = 1.0
because the performance difference among the solutions is
the largest. There is very slight difference in cell load levels
between the two cases. Numerically, the differences between
them are only 0.1% and 0.5% for average and maximum cell
load, respectively. This result is coherent with Figure 5(c) and
Figure 5(d). One can see that the patterns of the two pair
selection solutions are almost identical. Thus, pair filtering by
Lemma 1 is effective in reducing the number of candidate
pairs, with virtually no impact on performance.
In Figure 7 and Figure 8, we respectively evaluate PA(B-W)
and PA(B-SB) , combined with three power split schemes SP(Uni) ,
SP(FTPC) , and SP(Opt) . In SP(Opt) , we use the algorithm S PLIT to
compute the optimal power split for each pair. All of SP(Uni) ,
SP(FTPC) , and SP(Opt) are put into the framework of M-C ELL
but with fixed pair selection PA(B-W) or PA(B-SB) . In addition,
OMA(Opt) and NOMA(Opt) are also included for comparison as
baselines.
One can see that all the NOMA schemes outperform
OMA(Opt) . In Figure 7, with PA(B-W) , SP(FTPC) outperforms
SP(Uni) . SP(Opt) beats the other two. On one hand, there
is non-negligible gap in load performance between SP(Opt)
and NOMA(Opt) , even though in SP(Opt) , the power split is
optimal for the PA(B-W) pairing. Hence pair selection plays an
important role for NOMA performance. On the other hand,
SP(Opt) yields significantly load improvement compared to
OMA(Opt) , and we conclude that PA(B-W) is a good sub-optimal

pair selection for NOMA. Indeed, PA(B-W) pairs the UEs in
a greedy way, aiming at maximizing the diversity of channel
conditions of paired UEs. As shown in Figure 5(c), the optimal
pair selection has a similar trend. The difference is that optimal
pairing has a more “global view” than PA(B-W) . In Figure 8,
under PA(B-SB) , SP(Uni) , SP(FTPC) , and SP(Opt) improve the load
very slightly. All of the three are far from the global optimum
and the gap is large under high user demands. We conclude
that PA(B-SB) is not as effective as PA(B-W) in terms of network
load optimization.
As the overall conclusion, jointly optimizing power allocation and user pairing is important for the performance of
NOMA.
C. Convergence Analysis
We show the convergence performance of M-C ELL in Figure 9, for demands 0.3, 0.5, 0.7, and 1.0, respectively. Initially,
(0)
ρi = 1 (i ∈ C). We observe that M-C ELL converges very
fast. With higher demand, the convergence becomes slightly
faster. High accuracy is reached after a very few iterations. For
all the demands consider in the figure, even if we terminate
M-C ELL after a very few iterations, the obtained solution is
close to the optimum.
10−1
kρ(k) − ρ(k−1) k∞

power.
In summary, power allocation has considerably large influence on NOMA. Even if the UE pairs are optimally selected,
sub-optimal power allocations in NOMA have significant
deviation from optimal NOMA.

d = 1.0
d = 0.7
d = 0.4
d = 0.1

10−5
10−9
10−13
2

3

4

5

6

7

8

9

10

Iteration k in M-C ELL
Figure 9. This figure shows the norm k·k∞ in function of iteration k in
M-C ELL, under the uniform demands 0.1, 0.4, 0.7, and 1.0, respectively.

VIII. C ONCLUSIONS
This paper has investigated optimal resource management
in multi-cell NOMA, with power allocation and user pairing
being considered simultaneously. Joint optimization of both is
shown to be very important for NOMA performance. The proposed system model admits a mixed use of OMA and NOMA
for the users. Therefore, network architectures that support
various multiple access techniques can be analyzed under this
model. Finally, as for future work, the paper suggests that
mathematical tools in SIF are useful for analyzing multi-cell
NOMA. In summary, NOMA is a promising technique for
spectrum efficiency enhancement and cell load balancing.

13

A PPENDIX
We remark that the function (32) is first referred to in Section V-A and is used in the formulation (17). The two functions (33)
and (34) are first referred to in Section V-A and are used in the algorithm S PLIT.


w⊕ (ρ−i )e(c⊕u +c u ) + (w (ρ−i ) − w⊕ (ρ−i ))ec u
Cvu (c⊕u , c u , ρ−i ) = log
.
(32)
pi + w (ρ−i )

K
d⊕ /x
cK
u
⊕u xu > d⊕
h
i
d
/x
u
H⊕u (xu , ρ−i ) =
(33)
(pi +w (ρ−i ))/e ⊕
+w⊕ (ρ−i )−w (ρ−i )
log
Otherwise
w⊕ (ρ−i )

K K
d /x
hu
i c u xu > d
H u (xu , ρ−i ) =
(34)
pi +w (ρ−i )
Otherwise
log
d /xu
w⊕ (ρ−i )e

R EFERENCES
[1] L. You, L. Lei, D. Yuan, S. Sun, S. Chatzinotas, and B. Ottersten, “A
framework for optimizing multi-cell NOMA: Delivering demand with
less resource,” in 2017 IEEE GLOBECOM, 2017.
[2] S. M. R. Islam, N. Avazov, O. A. Dobre, and K. S. Kwak, “Powerdomain non-orthogonal multiple access (NOMA) in 5G systems: Potentials and challenges,” IEEE Communications Surveys Tutorials, vol. 19,
no. 2, pp. 721–742, 2017.
[3] W. Shin, M. Vaezi, B. Lee, D. J. Love, J. Lee, and H. V.
Poor, “Non-orthogonal multiple access in multi-cell networks: Theory,
performance, and practical challenges,” arXiv.org, 2016. [Online].
Available: https://arxiv.org/pdf/1611.01607.pdf
[4] Z. Ding, P. Fan, and H. V. Poor, “Impact of user pairing on 5G
nonorthogonal multiple-access downlink transmissions,” IEEE Transactions on Vehicular Technology, vol. 65, no. 8, pp. 6010–6023, 2016.
[5] “Evaluation methodologies for downlink multiuser superposition transmissions,” 3GPP, Tech. Rep. R1-153332, 2014.
[6] W. Shin, M. Vaezi, B. Lee, D. J. Love, J. Lee, and H. V. Poor,
“Coordinated beamforming for multi-cell MIMO-NOMA,” IEEE Communications Letters, vol. 21, no. 1, pp. 84–87, 2017.
[7] J. Kim, J. Koh, J. Kang, K. Lee, and J. Kang, “Design of user clustering
and precoding for downlink non-orthogonal multiple access (NOMA),”
in 2015 IEEE MILCOM, 2015, pp. 1170–1175.
[8] H. Tabassum, E. Hossain, and M. J. Hossain, “Modeling and analysis of
uplink non-orthogonal multiple access (NOMA) in large-scale cellular
networks using poisson cluster processes,” arXiv.org, 2016. [Online].
Available: http://arxiv.org/abs/1610.06995.pdf
[9] Y. Fu, Y. Chen, and C. W. Sung, “Distributed power control for the
downlink of multi-cell NOMA systems,” IEEE Transactions on Wireless
Communications, to appear.
[10] L. P. Qian, Y. Wu, H. Zhou, and X. Shen, “Joint uplink base station association and power control for small-cell networks with nonorthogonal multiple access,” IEEE Transactions on Wireless Communications, vol. 16, no. 9, pp. 5567–5582, 2017.
[11] B. Di, L. Song, Y. Li, and G. Y. Li, “Non-orthogonal multiple access
for high-reliable and low-latency V2X communications in 5G systems,”
IEEE Journal on Selected Areas in Communications, vol. 35, no. 10,
pp. 2383–2397, 2017.
[12] L. P. Qian, Y. Wu, H. Zhou, and X. Shen, “Dynamic cell association
for non-orthogonal multiple-access V2S networks,” IEEE Journal on
Selected Areas in Communications, vol. 35, no. 10, pp. 2342–2356,
2017.
[13] J. Zhu, J. Wang, Y. Huang, S. He, X. You, and L. Yang, “On optimal
power allocation for downlink non-orthogonal multiple access systems,”
IEEE Journal on Selected Areas in Communications, to appear.
[14] L. You and D. Yuan, “Joint CoMP-cell selection and resource allocation
in fronthaul-constrained C-RAN,” in 2017 WiOpt Workshop, 2017, pp.
1–6.
[15] L. Lei, D. Yuan, C. K. Ho, and S. Sun, “Optimal cell clustering and
activation for energy saving in load-coupled wireless networks,” IEEE
Transactions on Wireless Communications, vol. 14, no. 11, pp. 6150–
6163, 2015.
[16] I. Viering, M. Dottling, and A. Lobinger, “A mathematical perspective
of self-optimizing wireless networks,” in 2009 IEEE ICC, 2009, pp. 1–6.
[17] I. Siomina and D. Yuan, “Analysis of cell load coupling for LTE
network planning and optimization,” IEEE Transactions on Wireless
Communications, vol. 11, no. 6, pp. 2287–2297, 2012.

+w (ρ−i )−w⊕ (ρ−i )

[18] A. J. Fehske, I. Viering, J. Voigt, C. Sartori, S. Redana, and G. P.
Fettweis, “Small-cell self-organizing wireless networks,” Proceedings
of the IEEE, vol. 102, no. 3, pp. 334–350, 2014.
[19] L. You, D. Yuan, N. Pappas, and P. Värbrand, “Energy-aware wireless
relay selection in load-coupled OFDMA cellular networks,” IEEE Communications Letters, vol. 21, no. 1, pp. 144–147, 2017.
[20] L. You and D. Yuan, “Load optimization with user association in
cooperative and load-coupled LTE networks,” IEEE Transactions on
Wireless Communications, vol. 16, no. 5, pp. 3218–3231, 2017.
[21] I. Siomina, A. Furuskr, and G. Fodor, “A mathematical framework for
statistical QoS and capacity studies in OFDM networks,” in 2009 IEEE
PIMRC, 2009, pp. 2772–2776.
[22] K. Majewski and M. Koonert, “Conservative cell load approximation
for radio networks with Shannon channels and its application to LTE
network planning,” in 2010 Sixth Advanced International Conference on
Telecommunications, 2010, pp. 219–225.
[23] E. Pollakis, R. L. G. Cavalcante, and S. Staczak, “Base station selection for energy efficient network operation with the majorizationminimization algorithm,” in 2012 IEEE SPAWC, 2012, pp. 219–223.
[24] A. J. Fehske and G. P. Fettweis, “Aggregation of variables in load models
for interference-coupled cellular data networks,” in 2012 IEEE ICC,
2012, pp. 5102–5107.
[25] A. J. Fehske, H. Klessig, J. Voigt, and G. P. Fettweis, “Concurrent loadaware adjustment of user association and antenna tilts in self-organizing
radio networks,” IEEE Transactions on Vehicular Technology, vol. 62,
no. 5, pp. 1974–1988, 2013.
[26] C. K. Ho, D. Yuan, and S. Sun, “Data offloading in load coupled
networks: A utility maximization framework,” IEEE Transactions on
Wireless Communications, vol. 13, no. 4, pp. 1921–1931, April 2014.
[27] R. L. G. Cavalcante, S. Stanczak, M. Schubert, A. Eisenblaetter, and
U. Tuerke, “Toward energy-efficient 5g wireless communications technologies: Tools for decoupling the scaling of networks from the growth
of operating power,” IEEE Signal Processing Magazine, vol. 31, no. 6,
pp. 24–34, 2014.
[28] S. Tombaz, S. w. Han, K. W. Sung, and J. Zander, “Energy efficient
network deployment with cell DTX,” IEEE Communications Letters,
vol. 18, no. 6, pp. 977–980, 2014.
[29] B. Baszczyszyn, M. Jovanovic, and M. K. Karray, “Performance laws
of large heterogeneous cellular networks,” in 2015 WiOpt, 2015, pp.
597–604.
[30] C. K. Ho, D. Yuan, L. Lei, and S. Sun, “Power and load coupling
in cellular networks for energy optimization,” IEEE Transactions on
Wireless Communications, vol. 14, no. 1, pp. 509–519, 2015.
[31] R. L. G. Cavalcante, S. Staczak, J. Zhang, and H. Zhuang, “Low
complexity iterative algorithms for power estimation in ultra-dense load
coupled networks,” IEEE Transactions on Signal Processing, vol. 64,
no. 22, pp. 6058–6070, 2016.
[32] H. Klessig, D. hmann, A. J. Fehske, and G. P. Fettweis, “A performance
evaluation framework for interference-coupled cellular data networks,”
IEEE Transactions on Wireless Communications, vol. 15, no. 2, pp. 938–
950, 2016.
[33] R. L. G. Cavalcante, Y. Shen, and S. Staczak, “Elementary properties of
positive concave mappings with applications to network planning and
optimization,” IEEE Transactions on Signal Processing, vol. 64, no. 7,
pp. 1774–1783, 2016.
[34] Q. Liao, “Dynamic uplink/downlink resource management in flexible
duplex-enabled wireless networks,” in 2017 ICC Workshops, 2017, pp.
625–631.

14

[35] R. L. G. Cavalcante, M. Kasparick, and S. Staczak, “Max-min utility
optimization in load coupled interference networks,” IEEE Transactions
on Wireless Communications, vol. 16, no. 2, pp. 705–716, 2017.
[36] D. A. Awan, R. L. G. Cavalcante, and S. Stanczak, “A robust machine
learning method for cell-load approximation in wireless networks,”
arXiv.org, 2017. [Online]. Available: http://arxiv.org/abs/1710.09318
[37] D. Tse and P. Viswanath, Fundamentals of Wireless Communication.
Cambridge university press, 2005.
[38] G. Geraci, M. Wildemeersch, and T. Q. S. Quek, “Energy efficiency
of distributed signal processing in wireless networks: A cross-layer
analysis,” IEEE Transactions on Signal Processing, vol. 64, no. 4, pp.
1034–1047, 2016.
[39] M. Wildemeersch, T. Q. S. Quek, M. Kountouris, A. Rabbachin, and
C. H. Slump, “Successive interference cancellation in heterogeneous
networks,” IEEE Transactions on Communications, vol. 62, no. 12, pp.
4440–4453, 2014.
[40] R. D. Yates, “A framework for uplink power control in cellular radio
systems,” IEEE Journal on Selected Areas in Communications, vol. 13,
no. 7, pp. 1341–1347, 1995.
√
[41] S. Micali and V. V. Vazirani, “An O( V · |E|) algoithm for finding
maximum matching in general graphs,” in 21st Annual Symposium on
Foundations of Computer Science, 1980, pp. 17–27.
[42] H. R. Feyzmahdavian, M. Johansson, and T. Charalambous, “Contractive
interference functions and rates of convergence of distributed power
control laws,” IEEE Transactions on Wireless Communications, vol. 11,
no. 12, pp. 4494–4502, 2012.
[43] D. Huo, “Clarification on the wrap-around hexagon network structure,”
IEEE, Tech. Rep. C802.20-05/15, 2005.
[44] Y. Saito, A. Benjebbour, Y. Kishiyama, and T. Nakamura, “System-level
performance evaluation of downlink non-orthogonal multiple access
(NOMA),” in 2013 IEEE PIMRC, 2013, pp. 611–615.

