Joint Power Allocation and Beamforming for Energy-Efficient Two-Way Multi-Relay Communications
Zhichao Sheng Hoang D. Tuan Trung Q. Duong H. Vincent Poor





Joint Power Allocation and  Beamforming for Energy-Efficient  Two-Way Multi-Relay Communications
This work was supported in part  by the Australian Research Council’s Discovery Projects under Project DP130104617, in part by the U.K. Royal Academy of Engineering Research Fellowship under Grant RF14151422, and in part by
the U.S. National Science Foundation under Grants CCF-1420575 and CNS-1456793. 
    Zhichao Sheng, Hoang D. Tuan, Trung Q. Duong and H. Vincent  Poor
Zhichao Sheng and Hoang D. Tuan are with the Faculty of Engineering and Information Technology, University of Technology Sydney, Broadway, NSW 2007, Australia (email: kebon22@163.com, Tuan.Hoang@uts.edu.au)
Trung Q. Duong is with Queen's University Belfast, Belfast BT7 1NN, UK  (email: trung.q.duong@qub.ac.uk)
H. Vincent Poor is with the Department of Electrical Engineering, Princeton University, Princeton, NJ 08544, USA (e-mail: poor@princeton.edu)

    December 30, 2023
============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================





This paper considers the joint design of user power allocation and  relay beamforming in
relaying communications, in which multiple pairs of single-antenna users exchange information
with each other via multiple-antenna relays in two time slots. All users transmit their signals to the relays
in the first time slot while the relays broadcast the beamformed signals to all users in the second time slot.
The aim is to maximize the system's energy efficiency (EE) subject to
 quality-of-service (QoS) constraints in terms of exchange throughput requirements. The QoS constraints
are  nonconvex  with many nonlinear cross-terms, so finding a feasible point is already computationally challenging. The  sum throughput appears in the numerator while the total consumption power
appears in the denominator of the EE objective function. The former  is a  nonconcave function and the latter
is a nonconvex function, making fractional programming useless for EE optimization.
Nevertheless, efficient iterations of low complexity to obtain its optimized solutions are
developed. The performances of the multiple-user and multiple-relay networks under various scenarios are evaluated to show
the merit of the paper development.


Two-way relaying, information exchange (IE), energy efficiency (EE), quality-of-service (QoS),
relay beamforming, power allocation, joint optimization, path-following algorithms.





§ INTRODUCTION

Two-way relaying (TWR) <cit.> has been the focus of considerable research interest
in recent years due  to its potential in offering higher information exchange throughput for cognitive communications such as device-to-device (D2D) and machine-to-machine (M2M) communications <cit.>.  Unlike
the conventional one-way relaying, which needs four time slots for information exchange between a pair of users
(UEs), TWR needs just two time slots for this exchange <cit.>. In the first time slot, known as the multiple access (MAC) phase, all UEs simultaneously transmit their signals to the relays. In the second time slot, also known as the broadcast (BC) phase, the relays broadcast the beamformed signals to the all UEs.
Offering double fast communication, TWR obviously suffers from double
multi-channel interference as compared to one-way relaying  <cit.>. Both  optimal control of UEs' transmit power
and  TWR beamforming are thus very important in exploring the spectral efficiency of TWR.

There are various  scenarios of TWR considered in the literature. The most popular scenario is
single-antenna  relays serving a pair of single-antenna UEs <cit.>.
The typical problems are to design the  TWR  weights to either maximize the  throughput
or minimize the relay power subject to signal-to-interference-plus-noise ratio (SINR) constraints at UEs.
A branch-and-bound (BB) algorithm of global optimization <cit.> was used in <cit.>
for sum throughput maximization. Its computational complexity is already very high in very low-dimensional problems.
Semi-definite relaxation  of high computational complexity
followed by bisection  used in <cit.> works strictly under a single total relay power constraint.
Furthermore, the scenario of single-antenna relays serving multiple pairs of UEs was addressed in
<cit.> by  a polyblock algorithm of global optimization
<cit.> and in <cit.> by local linearization based iteration. The mentioned semi-definite relaxation
was also used in <cit.> in designing TWR beamformer for the scenario of
a multi-antenna relay serving multiple pairs of UEs. It should be realized that most of relay beamformer optimization
problems considered in all these works are not more difficult computationally than their one-way relaying counterparts,
which have been efficiently solved in <cit.>.

The fixed power allocation to UEs does not only miss the  opportunity of power distribution  within a network
but can also potentially increase interference to UEs of other networks <cit.>. The joint design of UE power allocation and TWR weights for single antenna relays serving a pair of single antenna UEs
to maximize the minimum SINR was considered in <cit.>.  Under the strict assumption  that
the complex channel gains from UEs to the relays are the same as those from the relays to the UEs,
its design is divided into two steps. The first step is to optimize the beamformer weights with UEs'
fixed power by sequential second-order convex cone programming (SOCP).   The second step
performs an exhaustive grid search for UE power allocation. Joint optimization of UE
precoding and relay beamforming for a multi-antenna relay serving a pair of multi-antenna UEs
could be successfully addressed only recently in <cit.>. An efficient computation for joint UE power
allocation and TWR beamforming to maximize the worst UE throughput  for multiple antenna relays
serving multiple pairs of single antenna UEs  was proposed
in <cit.>. The reader is also referred to <cit.> for joint UE and relay power allocation
in MIMO OFDM system of one relay and one pair of UEs.

Meanwhile, the aforementioned classical spectral efficiency (SE) in terms of high  throughput is now only
one among multiple driving forces for the development of the next generation communication networks (5G) <cit.>.
The energy consumption of communication systems has become sizable, raising environmental and economic
concerns <cit.>. Particularly, the network energy efficiency (EE) in  terms of
the ratio of the sum throughput and the total power consumption, which counts not only the transmission power but also the
drain efficiency of power amplifiers, circuit power and other power factors in supporting the network's activities,
is comprehensively  pushed forward in 5G to address these concerns <cit.>. EE in single-antenna TWR
has been considered in <cit.> for single-antenna OFDMA in assisting  multiple pairs of
single-antenna UEs and in <cit.> and <cit.> for multi-antenna relays in assisting a pair of multi-antenna UEs. Again, the
main tool for computational solution in these works is semi-definite relaxation, which not only
significantly increases the problem dimension but also performs unpredictably <cit.>. Also, the resultant
Dinkelbach's iteration of fractional programming invokes  a logarithmic function optimization, which is convex but still
computationally difficult with no available algorithm of polynomial complexity.


The above analysis of the state-of-the-art TWR motivates us to consider the joint  design of single-antenna UE power
allocation and TWR beamformers in a TWR network to maximize its EE subject to UE  QoS constraints.
We emphasize that both sum throughput maximization (for spectral efficiency) and EE maximization are meaningful
only in the context of UE QoS satisfaction, without which they will cause the UE service discrimination.
Unfortunately, to our best knowledge such UE QoS constraints were not addressed whenever they are nonconvex <cit.>. The nonconvexity of these QoS constraints implies that even finding their
feasible points is already computationally difficult. QoS constraints in terms of UEs' exchange throughput requirements
are much more difficult than that in terms of individual UE throughput requirements because the former cannot be
expressed in terms of individual SINR constraints  as the latter. To address the EE maximization problem, we first
develop a new computational method for UE exchange throughput requirement feasibility, which invokes only
simple convex quadratic optimization. A new path-following computational procedure for  computational solutions
of the EE maximization problem is then proposed.

The rest of the paper is organized as follows. Section
<ref> formulates  two optimization problems of EE maximization and UE QoS optimization.
Two  path-following algorithms are developed in Section <ref> and <ref> for their computation.
In contrast to fractional programming, these algorithms
invoke  only a simple convex quadratic optimization of low computational complexity at each iteration.
Section <ref> provides simulation results to verify the
performance of these algorithms. Finally, concluding remarks are
given in Section <ref>

Notation. Vectors and matrices are represented by boldfaced lowercase and uppercase, respectively.
 x(n) is the nth entry of vector x, while X(n,.) and X(n,m)  are the nth row and
  (n,m)th entry of a matrix X. ⟨x,
y⟩=x^Hy is the inner product between
vectors x and y. ||.|| is either the Euclidean vector squared norm or
the Frobenius matrix squared norm, and ℝ^N_+={x∈ℝ^N : x(n)≥ 0,n=1,2,…,N}.  X is the trace of matrix X
 and [X_m]_m=1,…,M is a block diagonal matrix with diagonal blocks X_m.
Lastly, x∼𝒞𝒩(x̅,R_x) means x is a vector of Gaussian random variables with mean x̅ and covariance R_x.





§ TWO-WAY RELAY NETWORKS WITH MULTIPLE MIMO RELAYS AND MULTIPLE SINGLE-ANTENNA USERS



Fig. <ref> illustrates a TWR network in which K pairs of single-antenna UEs exchange information with each other. Namely the kth UE (UE k) and the (K+k)th UE (UE K+k), with k=1,…,K, exchange information with
each other via
M relays designated as relay m, m=1,…,M,  equipped with N_R antennas.

Denote by s=(s_1,…,s_2K)^T ∈ℂ^2K  the vector of  information symbols transmitted by the UEs, whose
entries are independent and have unit energy, i.e., [ss^H]=I_2K, where 𝐈_2K is
the identity matrix of size (2K)× (2K).
Let h_ℓ,m∈ℂ^N_R be the vector of channels from UE ℓ to  relay m. The received signal at  relay m is

    r_m=∑_ℓ=1^2K√(p_ℓ)h_ℓ,m s_ℓ+n_R,m,

where  n_R,m∼𝒞𝒩(0,σ_R^2I_N_R) is the background noise, and p=(p_1,…,p_2K)^T ∈ℝ_+^2K represents the powers allocated to the UEs.

Relay m performs linear processing on the received signal by applying the beamforming matrix
W_m ∈^N_R × N_R. The beamformed signals are

    r_m,b=W_mr_m=∑_ℓ=1^2K√(p_ℓ)W_mh_ℓ,m s_ℓ+W_mn_R,m, m=1,…, M.

The transmit power at relay m is calculated as

    [||r_m,b||^2] = ∑_ℓ=1^2K p_ℓ ||W_mh_ℓ,m||^2+
       σ_R^2||W_m||^2.

Relay m transmits the beamformed signal r_m,b to the UEs.
Let g_m,k∈ℂ^N_R be the vector of channels from the relay m to UE k.
The received signal at  UE k is given by

    y_k   =   ∑_m=1^Mg_m,k^Tr_m,b+n_k
       =   ∑_m=1^Mg_m,k^T
    (∑_ℓ=1^2K√(p_ℓ)W_mh_ℓ,m s_ℓ+W_mn_R,m)+n_k,

where  n_k∼𝒞𝒩(0,σ_k^2) is the background noise, which can be written as

    [                                                                                                               y_k                                                                                                                 = √(p_χ(k))∑_m=1^Mg_m,k^T
    W_mh_χ(k),m s_χ(k)_ desired signal +
    √(p_k)∑_m=1^Mg_m,k^T
    W_mh_k,m s_k_ self-interference;                                                                                                                                    + ∑_m=1^Mg_m,k^T(∑_ℓ=1, ℓ≠ k, χ(k)^2K√(p_ℓ)W_mh_ℓ,m s_ℓ_ inter-pair interference+
    W_mn_R,m)+n_k. ]

In the above equation,  (k, χ(k)) is a pair of UEs that exchange information with each other, so

    χ(k)=[           K+k               k=1,…, K 
    k-K                 k=K+1,…, 2K ]

Assuming that the channel state information (CSI)  of the forward and backward channels  and the beamforming matrices
is available, UE k  effectively subtracts  the self-interference term in (<ref>) to have the SINR:

    γ_k(p, W)   =   p_χ(k)|∑_m=1^Mf_m,k^H
    W_mh_χ(k),m|^2/∑_ℓ=1, ℓ≠ k,χ(k)^2K p_ℓ|∑_m=1^Mf_m,k^H W_mh_ℓ,m|^2
    +σ_R^2∑_m=1^M||f_m,k^H W_m||^2+σ^2_k,

where f_m,k=g_m,k^*, so f_m,k^H=g_m,k^T. 

Under the definitions

    [                           _k,ℓ(W)≜∑_m=1^Mf_m,k^H W_mh_ℓ,m,; _k(W)≜[ f_1,k^H W_1 f_2,k^H W_2         ... f_M,k^H W_M ], ]

it follows that

    γ_k(p, W)   =    p_χ(k)|_k,χ(k)(W)|^2/
    [∑_ℓ=1, ℓ≠ k,χ(k)^2K p_ℓ|_k,ℓ(W)|^2+σ_R^2||_k(W)||^2 +σ^2_k].

In TWR networks, the UEs exchange information in bi-direction fashion in one time slot.
Thus, the  throughput at the kth UE pair
is defined by the following function of beamforming matrix W  and power allocation vector p:

    R_k(p, W)   =   ln ( 1+γ_k(p, W))+ln (1+γ_K+k(p, W)),     k=1,…,K.

Accordingly, the problem of maximizing the network EE subject to UE QoS constraints is
formulated as:


    max_W∈ℂ^N× N, p∈ℝ_+^2K   F_EE(𝐰,𝐩)≜∑_k=1^K[ln (1+ γ_k(p, W))+ln (1+γ_K+k(p, W))]/ζ(P^U_ sum(𝐩)+P^R_ sum(𝐩,
    𝐖))+MP^ R+2KP^ U
           0≤ p_k ≤ P_k^U,max,    k=1,…,2K, 
          ∑_k=1^2K p_k ≤ P_sum^U,max, 
          ∑_ℓ=1^2K p_ℓ||W_mh_ℓ,m||^2+
       σ_R^2||W_m||^2 ≤ P_m^A,max, m=1,...,M, 
          ∑_m=1^M(∑_ℓ=1^2K p_ℓ||W_mh_ℓ,m||^2+
       σ_R^2||W_m||^2)≤ P_sum^R,max, 
          R_k(p, W)≥ r_k, k=1,⋯, K,


where ζ, P^ R and P^ U are the reciprocal of drain efficiency of power amplifier,
the circuit powers of the relay and UE, respectively.  P^ R=N_R P_r and P_r is the circuit power for each antenna in relay.  (<ref>) is the exchange throughput QoS requirement for each pair of UEs.
Constraints (<ref>) and (<ref>) cap the transmit power of each UE k and each
relay m at predefined values P_k^U,max and P_m^A,max, respectively. On the other hand, constraints
(<ref>) and (<ref>) ensure that  the total transmit power of UEs and the total transmit
power of the relays not exceed the allowed power budgets P_sum^U,max and P_sum^R,max,
respectively.

Note  that (<ref>) is a very difficult nonconvex optimization problem because
the power constraints (<ref>) and (<ref>), the exchange throughput QoS
constraints (<ref>), and the objective function (<ref>)  are  nonconvex.
Moreover, the exchange throughput QoS constraints (<ref>) are much harder than the
typical individual throughput QoS constraints

    ln(1+γ_k(p,W))≥r̃_min, k=1,…, 2K,

which are equivalent to the computationally easier SINR constraints

    γ_k(p,W)≥ e^r̃_min-1, k=1,…, 2K.

Note that even finding  feasible point of the EE problem (<ref>) is already difficult as
it must be based on the following UEs QoS optimization problem


    max_W∈ℂ^N× N, p∈ℝ_+^2K φ(p,W)≜min_k=1,… ,K [ln (1+ γ_k(p, W))^1/r_k+ln (1+γ_K+k(p, W))^1/r_k]    
         (<ref>), (<ref>), (<ref>), (<ref>),


which is still highly nonconvex because
the objective function in (<ref>) is nonsmooth and nonconcave while the joint power constraints
(<ref>) and (<ref>)
in (<ref>) are nonconvex. Only a particular case of r_k≡ r_min was addressed in
<cit.>, under which  (<ref>) is then equivalent to the SINR multiplicative maximization


    max_W∈ℂ^N× N, p∈ℝ_+^2K   min_k=1,… ,K  [(1+γ_k(p, W))(1+γ_K+k(p, W)] 
          (<ref>), (<ref>), (<ref>), (<ref>),


which can be solve by d.c. (difference of two convex functions) iterations <cit.>.
To the authors' best knowledge there is no available computation for (<ref>) in general. The next
section is devoted to its solution.


§ MAXIMIN EXCHANGE THROUGHPUT OPTIMIZATION

To address (<ref>), our first major step is to transform the nonconvex constraints (<ref>) and (<ref>) to convex ones through variable change as follows.

Following <cit.>, make the variable change

    β_k=1/p_k^2≥ 0,    k=1,…,2K.

For α>0 and β>0, define the functions

    [                Ψ_k,ℓ(W, α, β)≜|_k,ℓ(W)|^2/√(αβ),    k,ℓ=1,…,2K,;                        Υ_k(W, α)≜||_k(W)||^2/√(α),    k=1,…,2K,; Φ_ℓ,m(W_m, α,β)≜||h_ℓ,m^HW_m||^2 /√(αβ),   ℓ=1,…,2K; m=1, …, M. ]

which are convex <cit.>.

The optimization problem (<ref>) can be equivalently rewritten as


    max_W∈ℂ^N× N, α∈ℝ_+^2K, β∈ℝ_+^2K f(W,α,
    β)≜min_k=1,… ,K 1/r_k[ln (1+|_k,K+k(W)|^2/√(α_kβ_K+k)).
       .+ln (1+|_K+k,k(W)|^2/√(α_K+kβ_k)) ]    
         ∑_ℓ=1, ℓ≠ k,χ(k)^2KΨ_k,ℓ(W, α_k, β_ℓ)
    +σ_R^2 Υ_k(W,α_k)+σ_k^21/√(α)_k≤ 1,  k=1,…,2K    
       β_k ≥1/(P_k^U,max)^2,    k=1,…,2K,    
        P_sum^U(β):=∑_k=1^2K1/√(β_k)≤ P_sum^U,max,    
       ∑_ℓ=1^2KΦ_ℓ,m(W_m,1,β_ℓ) +σ_R^2||W_m||^2  ≤ P_m^A,max,  m=1,…,M, 
       ∑_m=1^M[∑_ℓ=1^2KΦ_ℓ,m(W_m,1,β_ℓ) +σ_R^2||W_m||^2 ]
    ≤ P_sum^R,max.



One can see that γ_k(p, W) in (<ref>) is expressed in terms of functions in (<ref>) as

    γ_k(p, W)   =    |_k,χ(k)(W)|^2/1/p_χ(k)
    [∑_ℓ=1, ℓ≠ k,χ(k)^2KΨ_k,ℓ(W,1,1/p_ℓ^2)
    +σ_R^2Υ_k(W,1)+σ^2_k]

which is

    |_k,χ(k)(W)|^2/√(α_kβ_χ(k))

for

    β_k=1/p_k^2, k=1,…, 2K

and

    α_k=∑_ℓ=1, ℓ≠ k,χ(k)^2KΨ_k,ℓ(W,1,β_ℓ)
    +σ_R^2Υ_k(W,1)+σ^2_k.

Therefore,

    φ(𝐩,𝐖)=f(𝐖,α,β)

for α and β defined by (<ref>)-(<ref>), which  is also feasible for  (<ref>)
whenever  (𝐩,𝐖) is feasible for  (<ref>). We thus have proved that

    max (<ref>) ≤max (<ref>).

Note that (<ref>) is the same as

    ∑_ℓ=1, ℓ≠ k,χ(k)^2KΨ_k,ℓ(W,1,β_ℓ)
    +σ_R^2Υ_k(W,1)+σ^2_k≤√(α_k).

The point (𝐩,𝐖) with p_k=1/√(β)_k, k=1,…, 2K is feasible for (<ref>)
whenever (𝐖,α, β) is feasible for (<ref>). Using (<ref>) and (<ref>),
one can see

    f(W,α,β)≤φ(p,W),

implying

    max (<ref>) ≤max (<ref>).

The last inequality together with (<ref>) yield

    max (<ref>) ≤max (<ref>),

completing the proof of Theorem <ref>.

The benefit of expressing (<ref>) by (<ref>)  is that all constraints in the latter
are convex so the computational difficulty is concentrated in its objective function, which is lower bounded by
a concave function based on the following result.

 At (W^(κ),α^(κ),β^(κ)) it
is true that

    ln (1+|_k,χ(k)(W)|^2/√(α_kβ_χ(k)))   ≥   f^(κ)_k,χ(k)(W,α_k,β_χ(k))
       ≜   a_k,χ(k)^(κ)-b_k,χ(k)^(κ)√(α_k^(κ)β_χ(k)^(κ))
    [2{_k,χ(k)(W)(_k,χ(k)(W^(κ)))^*}
          -1/2|_k,χ(k)(W^(κ))|^2(α_k/α^(κ)_k
    +β_χ(k)/β^(κ)_χ(k))
    ]^-1

over the trust region

    2{_k,χ(k)(W)(_k,χ(k)(W^(κ)))^*}
    -1/2|_k,χ(k)(W^(κ))|^2(α_k/α^(κ)_k
    +β_χ(k)/β^(κ)_χ(k))>0,

for x_k,χ(k)^(κ)=|_k,χ(k)(W^(κ))|^2/√(α^(κ)_kβ^(κ)_χ(k)),

    [ a_k,χ(k)^(κ)=ln(1+x_k,χ(k)^(κ))+x_k,χ(k)^(κ)/x_k,χ(k)^(κ)+1>0,;                b_k,χ(k)^(κ)=(x_k,χ(k)^(κ))^2/x_k,χ(k)^(κ)+1>0. ]


We use the following inequalities with their proof given in the Appendix:

    ln(1+x) ≥ln(1+x̅) +x̅/x̅+1
    -x̅^2/x̅+11/x  ∀ x>0, x̅>0

and

    |x|^2/√(αβ)≥ 2{xx̅^*}/√(α̅β̅)
    -1/2|x̅|^2/√(α̅β̅)
    (α/α̅+β/β̅)  ∀ x∈ℂ, x̅∈ℂ,α>0,
    α̅>0, β>0,
    β̅>0,

which is the same as

    √(αβ)/|x|^2≤√(α̅β̅)/2{xx̅^*}
    -1/2|x̅|^2(α/α̅+β/β̅)  ∀ x∈ℂ, x̅∈ℂ,α>0,
    α̅>0, β>0,
    β̅>0.

Applying (<ref>) and (<ref>)
 for x=|_k,χ(k)(W)|^2/√(α_kβ_χ(k)) and x̅=|_k,χ(k)(W^(κ))|^2/√(α_k^(κ)β_χ(k)^(κ))
yields

    [                                             ln (1+|_k,χ(k)(W)|^2/√(α_kβ_χ(k)))                                                                              ≥                           a_k,χ(k)^(κ)-b_k,χ(k)^(κ)√(α_kβ_χ(k))/|_k,χ(k)(W)|^2;                                                                                                                                                             ≥ a_k,χ(k)^(κ)-b_k,χ(k)^(κ)√(α_k^(κ)β_χ(k)^(κ))
    [2{_k,χ(k)(W)(_k,χ(k)(W^(κ)))^*};                                                                                                  -1/2|_k,χ(k)(W^(κ))|^2(α_k/α^(κ)_k
    +β_χ(k)/β^(κ)_χ(k))
    ]^-1, ]

showing (<ref>). 

Accordingly, for a feasible (W^(κ), α^(κ), β^(κ)) of (<ref>)
found at the (κ-1)th iteration, the following convex optimization problem is solved
at the κth iteration to generate the next feasible  (W^(κ+1), α^(κ+1), β^(κ+1)):

    [ max_W∈ℂ^N× N, α∈ℝ_+^2K, β∈ℝ_+^2K
    f^(κ)(W,α,
    β)≜min_k=1,… ,K 1/r_k
    [f^(κ)_k,K+k(W,α_k,β_K+k)+f^(κ)_K+k,k(W,α_K+k,β_k)];                                                                                            (<ref>)-(<ref> ), (<ref>). ]

Note that

    f(W^(κ+1),α^(κ+1),β^(κ+1))≥ f^(κ)(W^(κ+1),α^(κ+1),
    β^(κ+1))

by (<ref>), and

    f(W^(κ),α^(κ),β^(κ))= f^(κ)(W^(κ),α^(κ),
    β^(κ))

because

    ln (1+|_k,χ(k)(W^(κ))|^2/√(α_k^(κ)β_χ(k)^(κ)))=f^(κ)_k,χ(k)(W^(κ),α_k^(κ),
    β_χ(k)^(κ)), k=1,…, 2K.

On the other hand, as (W^(κ),α^(κ),β^(κ)) and
(W^(κ+1),α^(κ+1),β^(κ+1)) are a feasible point and the optimal
solution of the convex optimization problem (<ref>), it is true that

    f^(κ)(W^(κ+1),α^(κ+1),
    β^(κ+1))>f^(κ)(W^(κ),α^(κ),
    β^(κ))

as far as (W^(κ),α^(κ),β^(κ))≠
(W^(κ+1),α^(κ+1),β^(κ+1)). The point (W^(κ+1),α^(κ+1),β^(κ+1)) is then better than (W^(κ),α^(κ),β^(κ)) because

    f(W^(κ+1),α^(κ+1),β^(κ+1))≥ f^(κ)(W^(κ+1),α^(κ+1)>f^(κ)(W^(κ),α^(κ),
    β^(κ))=f(W^(κ),α^(κ),β^(κ)).

Analogously to <cit.>, it can be shown that the sequence { (W^(κ),α^(κ),β^(κ)) } at least converges to a locally
optimal solution of the exchange throughput optimization problem (<ref>). The proposed Algorithm <ref>
for (<ref>) thus terminates after finitely many iteration, yielding an optimal solution
(𝐖^opt,α^opt,β^opt)
within tolerance ϵ>0. Then (𝐖^opt,p^opt) with p^opt=
(1/√(β^opt_1),…,1/√(β^opt_2K))^T
is accepted as the computational solution of the maximin
exchange throughput optimization problem (<ref>).


Before closing this section,  it is pointed out that the one-way
relay optimization in which UE k sends information to
UE  K+k  can be formulated as in (<ref>) by
setting γ_k=0 and  p_K+k=0 and thus can be directly solved by Algorithm <ref>.


§ ENERGY EFFICIENCY MAXIMIZATION

Return to consider the EE maximization problem (<ref>). It is worth noticing that the computational solution for the
QoS constrained sum throughput maximization problem

    max_W∈ℂ^N× N, p∈ℝ_+^2K∑_k=1^K[ln (1+ γ_k(p, W))+ln (1+γ_K+k(p, W))]
      (<ref>), (<ref>), (<ref>), (<ref>), (<ref>),

which is a particular case of (<ref>), is still largely open. As a by-product, our approach to computation for (<ref>)
is directly applicable to that for (<ref>).
Similarly to Theorem <ref>, it can be shown that (<ref>) is equivalently expressed by the following optimization
problem under the variable change (<ref>):


    max_W∈ℂ^N× N, α∈ℝ_+^2K, β∈ℝ_+^2K F(W,α,
    β)≜∑_k=1^K
    [ln (1+|_k,K+k(W)|^2/√(α_kβ_K+k))+ln (1+|_K+k,k(W)|^2/√(α_K+kβ_k))]/π(β,W)
       (<ref>)-(<ref>),
    R̃_k(𝐖,α,β)≜ln (1+|_k,K+k(W)|^2/√(α_kβ_K+k))+ln (1+|_K+k,k(W)|^2/√(α_K+kβ_k))≥ r_k, k=1,…, K,


where  the consumption power function π(β,W) is defined by

    π(β,W)≜∑_k=1^2Kζ/√(β_k)+ζ∑_m=1^M[∑_ℓ=1^2KΦ_ℓ,m(W_m,1,β_ℓ) +σ_R^2||W_m||^2 ]+MP^ R+2KP^ U.

In Dinkelbach's iteration based approach (see e.g. <cit.>), one aims to find through bisection
a value τ such that the optimal value of the following optimization problem is zero

    max_W∈ℂ^N× N, α∈ℝ_+^2K, β∈ℝ_+^2K[∑_k=1^K(ln (1+|_k,K+k(W)|^2/√(α_kβ_K+k))+ln (1+|_K+k,k(W)|^2/√(α_K+kβ_k))).
    .- π(β,W)]      (<ref>)-(<ref>). (<ref>).

Such τ obviously is the optimal value of (<ref>). However, for each τ,  (<ref>) is still
nonconvex and as computationally difficult as the original optimization problem (<ref>).
There is no benefit to use (<ref>).

To address computation for (<ref>) involving the  nonconcave objective function F(W,α,
β) and the nonconvex constraint (<ref>),
we  will explore the following inequality for positive quantities, whose proof is given
in the Appendix:

    ln(1+x)/t   ≥    2ln(1+x̅)/t̅+x̅/t̅(x̅+1)-x̅^2/(x̅+1)t̅1/x-
    ln(1+x̅)/t̅^2t 
            ∀ x>0, x̅>0, t>0, t̅>0.

The right-hand-side (RHS) of (<ref>) is a concave function on the interior domain of R^2_+
and agrees with the left-hand-side (LHS) at (x̅,t̅).

Suppose that (W^(κ),α^(κ),β^(κ)) is a feasible point of (<ref>)
found from the (κ-1)th iteration. Applying (<ref>) for

    x=|_k,χ(k)(W)|^2/√(α_kβ_χ(k)), t=π(β,W)

and

    x̅=|_k,χ(k)(W^(κ))|^2/√(α_k^(κ)β_χ(k)^(κ)),
    t̅=π(β^(κ),W^(κ))

yields

    ln (1+|_k,χ(k)(W)|^2/√(α_kβ_χ(k)))/π(β,W)   ≥   F_k,χ(k)^(κ)(W,α_k,
    β)
       ≜   p_k,χ(k)^(κ)- q_k,χ(k)^(κ)√(α^(κ)_kβ^(κ)_χ(k))[2{_k,χ(k)(W)(_k,χ(k)(W^(κ)))^*}
          -1/2|_k,χ(k)(W^(κ))|^2(α_k/α^(κ)_k
    +β_χ(k)/β^(κ)_χ(k))
    ]^-1
          -r_k,χ(k)^(κ)π(β,W),

over the trust region (<ref>),  for x_k,χ(k)^(κ)=|_k,χ(k)(W^(κ))|^2/√(α^(κ)_kβ^(κ)_χ(k)),

    [                                                           t^(κ)                                                               =                                                 π(β^(κ),W^(κ)),;                                                    p_k,χ(k)^(κ)                                                               = 2ln(1+x_k,χ(k)^(κ))/t^(κ)+
    x_k,χ(k)^(κ)/t^(κ)(x_k,χ(k)^(κ)+1)>0;                                                    q_k,χ(k)^(κ)                                                               =                       (x_k,χ(k)^(κ))^2/(x_k,χ(k)^(κ)+1)t^(κ)>0,;                                                    r_k,χ(k)^(κ)                                                               =                                 ln(1+x_k,χ(k)^(κ))/(t^(κ))^2>0. ]

Accordingly, the following convex optimization problem is solved at the κth iteration  to generate the next
iterative point (W^(κ+1),α^(κ+1),β^(κ+1)):


    max_W∈ℂ^N× N, α∈ℝ_+^2K, β∈ℝ_+^2K F(W,α,
    β)≜∑_k=1^K
    [F^(κ)_k,K+k(W,α_k,β)+F^(κ)_K+k,k(W,α_K+k,β)]
       (<ref>)-(<ref> ), (<ref>), (<ref>),
    
    f^(κ)_k,K+k(W,α_k,β_K+k) +f^(κ)_K+k,k(W,α_K+k,β_k)
    ≥ r_k, k=1,…, K,


where f^(κ)_k,χ(k) are defined from (<ref>). Since,

    f^(κ)_k,K+k(W,α_k,β_K+k) +f^(κ)_K+k,k(W,α_K+k,β_k)≤R̃_k(𝐖,α,β)

it follows that the feasibility of the nonconvex constraint (<ref>) in (<ref>) is guaranteed by that
of (<ref>). Also,

    f^(κ)_k,K+k(W^(κ),α_k^(κ),β_K+k^(κ)) +f^(κ)_K+k,k(W^(κ),α_K+k^(κ),β_k^(κ))=
    R̃_k((W^(κ),α^(κ),β^(κ)))≥ r_min

because (W^(κ),α^(κ),β^(κ)) is feasible for (<ref>) and thus
feasible for (<ref>). Therefore, the convex optimization problem (<ref>) is always feasible. Analogously
to the previous section, the sequence {(W^(κ),α^(κ),β^(κ))} is seen
convergent at least to a locally optimal solution of problem (<ref>) and as thus the proposed
Algorithm  terminates after finitely many iteration, yielding an optimal solution
(W^(opt),α^opt,β^opt) within tolerance ϵ. Then
(W^opt,p^opt) with p^opt=(1/√(β_1^opt),...,1/√(β_2K^opt))^T is accepted as
the computational solution of the EE maximization problem (<ref>).


To find an initial feasible point (W^(0), α^(0), β^(0)) for Algorithm <ref>,
we use Algorithm <ref> for
computing (<ref>), which terminates upon

    [min_k=1,… ,K R_k(W^(κ), α^(κ), β^(κ))/r_k]≥ r_min

to satisfy (<ref>).

For the QoS constraints (<ref>) instead of (<ref>), which by the variable change (<ref>)
are equivalent to the following constraints

    |_k,χ(k)(W)|^2-(e^r̃_min-1)√(α_kβ_χ(k))≥ 0, k=1,…, 2K.

The LHS of (<ref>) is a convex functions, so (<ref>) is called reverse convex <cit.>, which can
be easily innerly approximated by linear approximation of the LHS at (W^(κ),α_k^(κ),
β_χ(k)^(κ)).

Remark. To compare the energy efficiency with one-way communication we need to revisit the one-way model: the users
{1,⋯,K} send symbols (s_1,⋯,s_K)^T∈ℂ^K via the relays in the first stage and the
users {K+1,...,2K} send symbols (s_K+1,⋯,s_2K)^T∈ℂ^K via the relays in the second state.
Denote by W^1_m and W^2_m the beamforming matrices for the received signals from the users
{1,⋯,K} and {K+1,⋯,2K}, respectively.  The transmit power at relay m in forwarding signals to users {K+1,⋯, 2K} in the first stage is

    P_m^A,1(p^1,W^1_m) = ∑_ℓ=1^K p_ℓ ||W^1_mh_ℓ,m||^2+
       σ_R^2||W^1_m||^2,

and the transmit power at relay m in forwarding signals to users {1,⋯, K} in the second stage is

    [ P_m^A,2(p^2,W^2_m) = ∑_ℓ=1^K p_ℓ+K ||W^2_mh_ℓ,m||^2+
       σ_R^2||W^2_m||^2. ]

Therefore, the power constraint  at  relay m is

    ∑_ℓ=1^K
     ( p_ℓ||W^1_mh_ℓ,m||^2+ p_ℓ+K||W^2_mh_ℓ+K,m||^2)+
       σ_R^2(||W^1_m||^2+||W^2_m||^2)≤ P_m^A,max, 
    
       m=1,...,M.

The total  power constraint is

    ∑_m=1^M[∑_ℓ=1^K
     ( p_ℓ||W^1_mh_ℓ,m||^2+ p_ℓ+K||W^2_mh_ℓ+K,m||^2)+
       σ_R^2(||W^1_m||^2+||W^2_m||^2)
    ]≤  P_sum^R,max.

Accordingly, for k=1,⋯, K, the SINR at UEs can be calculated as:

    γ̃_K+k(p^1, W^1)   =   p_k|∑_m=1^Mf_m,K+k^H
    W^1_mh_k,m|^2/∑_ℓ=1, ℓ≠ k^K p_ℓ|∑_m=1^Mf_m,K+k^H W^1_mh_ℓ,m|^2
    +σ_R^2∑_m=1^M||f_m,K+k^H W^1_m||^2+σ^2_K+k
       =    p_k|_K+k,k(W^1)|^2/
    [∑_ℓ=1, ℓ≠ k^K p_ℓ|_K+k,ℓ(W^1)|^2
          +σ_R^2||_K+k(W^1)||^2 +σ^2_K+k]

and

    γ̃_k(p^2, W^2)   =   p_K+k|∑_m=1^Mf_m,k^H
    W^2_mh_K+k,m|^2/∑_ℓ=1, ℓ≠ k^K p_K+ℓ|∑_m=1^Mf_m,k^H W^2_mh_K+ℓ,m|^2
    +σ_R^2∑_m=1^M||f_m,k^H W^2_m||^2+σ^2_k
       =    p_K+k|_k,K+k(W^2)|^2/
    [∑_ℓ=1, ℓ≠ k^K p_K+ℓ|_k,K+ℓ(W^2)|^2
          +σ_R^2||_k(W^2)||^2 +σ^2_k].

The EE maximization problem is then formulated as


    max_W^1∈ℂ^N× N, W^2∈ℂ^N× N, p∈ℝ_+^2K   1/2∑_k=1^K[ln (1+ γ̃_k(p^2, W^2))+ln (1+γ̃_K+k(p^1, W^1))]/ζ(P^U_ sum(𝐩)+P̃^R_ sum(𝐩,
    𝐖))+2MP^ R+2KP^ U
          (<ref>), (<ref>), (<ref>), (<ref>), 
          ln (1+ γ̃_k(p^2, W^2))+ln (1+γ̃_K+k(p^1, W^1))
    ≥ r_min, 
          k=1,⋯, K.


The pre-log factor 1/2 in the numerator of (<ref>) is to account for two stages needed in
communicating s_1,⋯, 2K and the non transmission power consumption at the relays  2MP^ R to
reflect the fact that the relays have to transmit twice.



§ NUMERICAL RESULTS

This sections evaluates the proposed algorithms  through the simulation.
The channels in the receive signal equations (<ref>) and (<ref>) are assumed Rayleigh fading,
which are modelled by independent circularly-symmetric complex Gaussian random variables
with zero mean and unit variance, while  the background noises n_R,m
and n_k  are also normalized, i.e., σ_R^2=σ_k^2=1. The computation tolerance
for terminating the Algorithm is ϵ=10^-4.
The numerical results are averaged over 1,000 random channel realizations.

Without loss of generality, simply set P_k^U, max≡ P^U, max, and P_n^A, max≡ P^A, max, P_sum^U,max=KP^U, max, and P_sum^R,max=MP^A, max/2.
P^U, max is fixed at 10 dBW  but the relay power budget P_sum^R,max
varies from 0 to 30 dBW.    The drain efficiency of power amplifier 1/ζ is set
to be 40 %. As in <cit.>, the circuit powers of each antenna in relay and UE are 0.97 dBW and -13 dBW, respectively. We consider the scenarios of K ∈{1,2,3} pairs
and (M,N_R) ∈{(1, 8), (2, 4), (4, 2)}, i.e. the total number of antennas is fixed at 8 but the number of relays
is M∈{1, 2, 4}.



 §.§ Maximin exchange throughput optimization

This subsection analyses  the exchange throughput achieved by  TWR.
The jointly optimal power and relay beamforming is referred  as OP-OW, while
the optimal beamforming weights with UE equal power allocation are referred as OW. The initial points for Algorithm
<ref> is chosen from the OW solutions.
To compare with the numerical result of <cit.>, r_k≡ 1 is set for (<ref>).

Fig. <ref>, <ref> and <ref> plot the
achievable  minimum pair exchange throughput versus the  relay power budget P_sum^R,max with K ∈{1,2,3}.
The improvement by OPOW over OW  is significant for K=2 and K=3.
The throughput gain is more considerable  with higher P_sum^R,max. It is also observed that
using less relays achieves better throughput.
 The result in Fig. <ref> for K=2 is consistent with that in <cit.>.






Table <ref> provides the average number of iterations of Algorithm <ref>. As can be observed,
Algorithm <ref>  converges in less than 24 iterations in all considered scenarios.




 §.§ EE maximization

This subsection examines the performance of energy efficiency achieved  by Algorithm <ref>.
r_k in (<ref>) is set as the half of the optimal value obtained by Algorithm <ref>.
Firstly, the simulation results presented in Fig. <ref>, <ref> and <ref> are for
K=2 and (M,N_R) ∈{(1, 8), (2, 4), (4, 2)}.
Fig. <ref> compares the EE performance achieved by TWR, one-way relaying
and TWR with UE fixed equal power allocation, which is labelled as "other method".
It is clear from Fig. <ref> that TWR clearly and significantly outperforms one-way relaying and
TWR UE equal power allocation.

Under small transmit power regime, the power consumption in the denominator is dominated by the circuit power and
the EE is maximized by maximizing the sum throughput in the numerator. As such all the EE, the sum
throughput and transmit power increase in Fig. <ref>, <ref> and
<ref> in the relay power budget P_sum^R,max.

However, under larger transmit power regime, where the denominator of EE is dominated by the actual transmit power,
the EE becomes to be  maximized by minimizing the transmit power in the denominator,
which saturates when beyond a threshold. When the transmit power saturates in Fig. <ref>, both the
EE and the sum throughput accordingly saturate in Fig. <ref> and <ref>.
It is also observed that for a given  relay power budget and a given number of total antennas in all relays, the configuration
with less relays is superior to the ones with more relays. This is quite expected since the configuration
with less relays achieves higher throughput than the ones with more relays.
Table <ref> shows that Algorithm <ref> converges
in less than 26 iterations.

Similar comparisons are provided in Fig. <ref> and <ref> for K=1 and K=3, respectively
with the superior  EE performances of TWR over one-way relaying and TWR UE equal power allocation observed.




















§ CONCLUSIONS

Joint UE power allocation and relay beamforming  to either satisfy a UE's QoS requirement or
maximize the energy efficiency of TWR serving multiple UEs is a very difficult nonconvex optimization problem.
This paper has developed two path-following computational procedures for their solutions, which invoke a simple
convex quadratic program of low computational complexity at each iteration.
Simulation results have confirmed their rapid convergence. We have shown that TWR achieves  much higher energy-efficiency
than its one-way relaying counterpart in all considered scenarios.


§ APPENDIX: PROOF FOR INEQUALITIES (<REF>), (<REF>) AND (<REF>)

The function ψ_1(z)=ln (1+z^-1) is convex on the domain z>0 <cit.>. Therefore <cit.>

    ψ_1(z)≥ψ(z̅) +∇ψ_1(z̅)(z-z̅)   ∀ z>0, z̅>0,

which is seen as

    ln (1+z^-1)≥ln(1+z̅^-1)+1/z̅+1-z/(z̅+1)z̅  ∀ z>0, z̅>0.

The inequality (<ref>) is obtained by substituting x=z^-1 and x̅=z̅^-1 into (<ref>).

Next, the function ψ_2(x,α,β)=|x|^2/√(αβ) is convex on the domain
z∈ℂ, α>0 and β>0 <cit.>, so again

    ψ_2(x,α,β)≥ψ_2(x̅,α̅,β̅)+∇ψ_2(x̅,α̅,β̅), (x,α,β)-(x̅,α̅,β̅),

which is seen as (<ref>).

Finally, by checking its Hessian, the function ψ_3(z,t)=(ln (1+z^-1))/t is seen to be
convex on the interior of R^2_+.
Therefore,

    ψ_3(z,t)≥ψ_2(z̅,t̅)+∇ψ_3(z̅,t̅), (z,t)-(z̅,t̅)  ∀ z>0, z̅>0, t>0, t̅>0,

which is seen as


    ln(1+z^-1)/t≥ 2ln(1+z̅^-1)/t̅+1/(z̅+1)t̅-
    z/(z̅+1)z̅t̅-ln(1+z̅)/t̅^2t.

The inequality (<ref>) follows by substituting x=z^-1 and x̅=z̅^-1 into the last inequality.
10

Amarasuriya12
G. Amarasuriya, C. Tellambura, and M. Ardakani, “Two-way amplify-and-forward
  multiple-input multiple-output relay networks with antenna selection,” IEEE J. Sel. Areas Commun., vol. 30, pp. 1513–1529, Sep. 2012.

Chung12
H. Chung, N. Lee, B. Shim, and T. Oh, “On the beamforming design for MIMO
  multipair two-way relay channels,” IEEE Trans. Vehicular Technology,
  vol. 16, pp. 3301–3306, Sep. 2012.

Getal13
D. Gunduz, A. Yener, A. Goldsmith, and H. V. Poor, “The multiway relay
  channel,” IEEE Trans. Inf. Theory, vol. 59, pp. 51–63, Jan 2013.

AWM14
A. Asadi, Q. Wang, and V. Mancuso, “A survey device-to-device communication in
  cellular networks,” IEEE Commun. Surveys & Tut., vol. 16,
  pp. 1801–1819, 4th Qart. 2014.

Liuetal15
J. Liu, N. Kato, J. Ma, and N. Kadowaki, “Device-to-device communication in
  LTE-advanced networks: A survey,” IEEE Commun. Surveys Tuts.,
  vol. 17, pp. 1923–1940, 4th Quart. 2015.

ZhangL09
R. Zhang, Y.-C. Liang, C. C. Chai, and S. Cui, “Optimal beamforming for
  two-way multi-antenna relay channel with analogue network coding,” IEEE
  J. Sel. Areas Commun., vol. 27, pp. 699–712, Jul. 2009.

Zeng11
M. Zeng, R. Zhang, and S. Cui, “On design of collaborative beamforming for
  two-way relay networks,” IEEE Trans. Signal Prosess., vol. 59,
  pp. 2284–2295, May 2011.

Chalise10
B. K. Chalise and L. Vandendorpe, “Optimization of MIMO relays for
  multipoint-to-multipoint communications: Nonrobust and robust designs,” IEEE Trans. Signal Prosess., vol. 58, pp. 6355–6368, Dec. 2010.

Phan13
A. H. Phan, H. D. Tuan, H. H. Kha, and H. H. Nguyen, “Iterative D.C.
  optimization of precoding in wireless MIMO relaying,” IEEE Trans.
  Wireless Commun., vol. 12, pp. 1617–1627, Apr. 2013.

Cheng11
W. Cheng, M. Ghogho, Q. Huang, D. Ma, and J. Wei, “Maximizing the sum-rate of
  amplify-and-forward two-way relaying networks,” IEEE Signal Prosess.
  Lett., vol. 18, pp. 635–638, Nov. 2011.

Amirani12
M. Zaeri-Amirani, S. Shahbazpanahi, T. Mirfakhraie, and K. Ozdemir,
  “Performance tradeoffs in amplify-and-forward bidirectional network
  beamforming,” IEEE Trans. Signal Prosess., vol. 60, pp. 4196–7209,
  Aug. 2012.

Tuybook
H. Tuy, Convex Analysis and Global Optimization (second edition).
Springer, 2016.

Zhang12
J. Zhang, F. Roemer, M. Haardt, A. Khabbazibasmenj, and S. A. Vorobyov, “Sum
  rate maximixation for multi-pair two-way relaying single-antenna amplify and
  forward relays,” in Proc. IEEE Int. Conf. Acoust Speech Signal Prosess.
  (ICASSP), Mar. 2012.

Tao12
M. Tao and R. Wang, “Linear precoding for multi-pair two-way MIMO relay
  systems with max-min fairness,” IEEE Trans. Signal Prosess., vol. 60,
  pp. 5361–5370, Oct. 2012.

Wang12
R. Wang, M. Tao, and Y. Huang, “Linear precoding designs for
  amplify-and-forward multiuser two-way relay systems,” IEEE Trans.
  Wireless Commun., vol. 11, pp. 4457–4469, Dec. 2012.

KRVH12
A. Khabbazibasmenj, F. Roemer, S. A. Vorobyov, and M. Haardt, “Sum-rate
  maximization in two-way AF MIMO relaying: Polynomial time solutions to a
  class of DC programming problems,” IEEE Trans. Signal Prosess.,
  vol. 60, pp. 5478–5493, Dec. 2012.

Khandaker11
M. R. A. Khandaker and Y. Rong, “Joint power control and beamforming for
  peer-to-peer MIMO relay systems,” in Proc. IEEE Int. Conf. Wireless
  Commun. Signal Prosess. (WCSP), Nov. 2011.

ChengPesavento12
Y. Cheng and M. Pesavento, “Joint optimization of source power allocation and
  distributed relay beamforming in multiuser peer-to-peer relay networks,”
  IEEE Trans. Signal Prosess., vol. 60, pp. 2962–2973, Jun. 2012.

Jing12
Y. Jing and S. Shahbazpanahi, “Max-min optimal joint power control and
  distributed beamforming for two-way relay networks under per-node power
  constraints,” IEEE Trans. Signal Prosess., vol. 60, pp. 6576–6589,
  Dec. 2012.

Raetal14
U. Rashid, H. D. Tuan, H. H. Kha, and H. H. Nguyen, “Joint optimization of
  source precoding and relay beamforming in wireless MIMO relay networks,”
  IEEE Trans. Commun., vol. 62, no. 2, pp. 488–499, 2014.

Khetal13
H. H. Kha, H. D. Tuan, H. H. Nguyen, and H. H. M. Tam, “Joint design of user
  power allocation and relay beamforming in two-way mimo relay networks,” in
  Proc. 7th Inter. Conf. Signal Process. Commun. Syst. (ICSPCS),
  pp. 1–6, 2013.

TNT17
H. D. Tuan, D. T. Ngo, and H. H. M. Tam, “Joint power allocation for
  MIMO-OFDM full-duplex relaying communications,” EURASIP J. Wireless
  Commun. Networking, 2017, DOI 10.1186/s13638-016-0800-4.

BJDO14
E. Bjornson, E. A. Jorswieck, M. Debbah, and B. Ottersten, “Multiobjective
  signal processing optimization: The way to balance conflicting metrics in
  5G systems,” IEEE Signal Process. Mag., vol. 31, pp. 14–23, Nov.
  2014.

Feetal11
A. Fehske, G. Fettweis, J. Malmodin, and G. Biczok, “The global footprint of
  mobile communications: The ecological and economic perspective,” IEEE
  Commun. Mag., vol. 49, pp. 55–72, Aug. 2011.

Buetal16
S. Buzzi, C.-L. I, T. E. Klein, H. V. Poor, C. Yang, and A. Zappone, “A survey
  of energy-efficient techniques for 5G networks and challenges ahead,” IEEE J. Select. Areas Commun., vol. 34, pp. 697–709, Apr. 2016.

XLL15
C. Xiong, L. Lu, and G. Y. Li, “Energy-efficient OFDMA-based two-way
  relay,” IEEE Trans. Commun., vol. 63, pp. 3157–3169, Sep. 2015.

ZH15
J. Zhang and M. Haardt, “Energy efficient two-way non-regenerative relaying
  for relays with multiple antennas,” IEEE Signal Process. Lett.,
  vol. 22, pp. 1079–1083, Aug. 2015.

WLTP16
Z. Wang, L. Li, H. Tian, and A. Paulraj, “User-centric precoding designs for
  the non-regenerative MIMO two-way relay systems,” IEEE Commun.
  Lett., vol. 20, pp. 1935–1938, Oct. 2016.

Phan12
A. H. Phan, H. D. Tuan, H. H. Kha, and D. T. Ngo, “Nonsmooth optimization for
  efficient beamforming in cognitive radio multicast transmission,” IEEE
  Trans. Signal Process., vol. 60, pp. 2941–2951, Jun. 2012.

KTN12
H. H. Kha, H. D. Tuan, and H. H. Nguyen, “Fast global optimal power allocation
  in wireless networks by local D.C. programming,” IEEE Trans. on
  Wireless Commun., vol. 1, pp. 510–515, Feb. 2012.

DM08
B. Dacorogna and P. Marechal, “The role of perspective functions in convexity,
  polyconvexity, rank-one convexity and separate convexity,” J. of Convex
  Analysis, vol. 15, no. 2, pp. 271–284, 2008.

TTN16
H. H. M. Tam, H. D. Tuan, and D. T. Ngo, “Successive convex quadratic
  programming for quality-of-service management in full-duplex MU-MIMO
  multicell networks,” IEEE Trans. Comm., vol. 64, pp. 2340–2353, Jun.
  2016.

Xiong2012
C. Xiong, G. Y. Li, S. Zhang, Y. Chen, and S. Xu, “Energy-efficient resource
  allocation in OFDMA networks,” IEEE Trans. Commun., vol. 60,
  pp. 3767–3778, Dec. 2012.

Taetal16
H. M. T. Ho, H. D. Tuan, D. Ngo, T. Q. Duong, and V. Poor, “Joint load
  balancing and interference management for small-cell heterogeneous networks
  with limited backhaul capacity,” IEEE Trans. Wireless Commun.,
  vol. PP, no. 99, pp. 1–1, 2016.



