1 Introduction

Recently, unmanned aerial vehicles (UAVs) have received dramatical attentions in wireless communication networks since they have the capability of flexible mobility [14]. Compared to the traditional wireless terminals, UAVs can be regarded as mobile relays or base stations. Therefore, UAV plays a critical role in the performance improvement and coverage extensiveness of the wireless communication network [5, 6]. In particular, UAV has huge potential in the intractable transmission environment, such as mountainous and emergency zones. Hence, the UAV which utilized as a relay has been widely investigated in the wireless network system. Specially, in [7], Song et al. investigated the instantaneous rate maximization problem for UAV-enabled decode-and-forward (DF) relay system by jointly optimizing the precoding matrix and power allocation. Ji et al. [8] studied the secrecy performance analysis in a UAV DF relay system with energy harvesting. In [9], HU et al. considered the energy consumption minimization problem in UAV-assisted relay system, where UAV is regarded as a computation server for user equipment (UE) and forwards the result to the access point. Compared with one-way relay system, two-way relay system has higher spectrum efficiency. In [10], Li et al. provided a UAV-enable two-way relay system, where a set of UEs was considered. The UAV trajectory and the transmit power of all terminals were jointly optimized to maximize the sum rate of UEs.

In addition to the service coverage, the energy constraints of terminals also have a critical effect on the effectiveness and reliability of wireless networks [1114]. In order to tackle this problem, simultaneous wireless information and power transfer (SWIPT) which can take advantage of the broadcast character of radio-frequency (RF) signal is adopted to enhance the performance of wireless communication networks [1519]. In [20], two classical and feasible schemes, i.e., time switching (TS) and power splitting (PS) are proposed for SWIPT technique. The TS-based design switches over time between energy harvesting (EH) and information decoding (ID) processing, while the PS-based one splits the received RF signals into two power streams, one for ID and the other for EH. A multiple-input multiple-output (MIMO) relaying system with TS protocol was considered in [21], where the rate was maximized by designing the robust beamforming matrixes of transceivers under energy constraints. In [22], Wang et al. investigated a multi-antenna relay communication system with a TS-powered relay node. Furthermore, a robust precoding scheme that can maximize the sum-rate of the multiuser relay system has been proposed in [22]. Subsequently, in [23], a joint design of PS factor and the beamforming matrixes at the transceivers for MIMO two-way relay system based on maximizing the energy efficiency was introduced.

Inspired by the advantage of SWIPT, more and more researchers and scientists are attempting to adopt the SWIPT technique to UAV-enabled relay system for improving the system performance [2426]. In [24], Wang et al. provided a UAV-enable non-orthogonal multiple access (NOMA) relay system with SWIPT, where the nonlinear energy harvesting model was considered. The PS factor and beamforming vectors were jointly optimized to maximize the rate of UEs while still guaranteeing security. Yin et al. [25] considered the throughput maximization problem in a UAV relay system with time-sharing mechanism, where the UAV relay is capable of SWIPT. In [26], the authors investigated a millimeter-wave UAV relay system with SWIPT to improve the secrecy performance, and the security rate of NOMA and orthogonal multiple access (OMA) were derived. However, the robust beamforming and SWIPT design for UAV-enabled relay system, where the UEs are capable of PS scheme, has not been well addressed.

2 Methods

In order to address the above practical issue, we investigate the UAV-enabled two-way relaying system with SWIPT. The main contributions are fourfold.

  • We propose a UAV-enabled two-way relaying system with SWIPT. The system model distinguishes from the existing relaying systems, since the UEs are capable of PS protocol and multi-antennas are equipped at UEs.

  • A novel optimization problem is formulated by maximizing the sum remained energy of two UEs subject to both UAV and UEs transmission power, while guaranteeing the minimum signal-to-noise ratio (SNR) requirement at UEs.

  • We propose an alternating optimization (AO) scheme and a low-complexity scheme based on the generalized singular value decomposition (GSVD) for the formulated optimization problem.

  • Simulation and numerical results are conducted to evaluate the robustness and effectiveness of two proposed schemes. Besides, the location ratio has been investigated to analyze the performance of the relay system.

The rest of the paper is outlined as follows. In Section 3, the UAV-enabled two-way relaying system with two PS UEs is described in detail. Section 4 concentrates on the design of the proposed AO scheme for maximizing the sum remained energy of two UEs. Section 5 provides the low-complexity scheme based on GSVD. In Section 6, numerical results are conducted to validate the performance of two proposed schemes. Finally, Section 7 gives the conclusion.

3 System model and problem formulation

3.1 System model

As shown in Fig. 1, we consider a UAV-enable two-way relay system where a UAV relay equipped with N antennas helps two UEs to exchange information and two UEs, i.e., UE 1 and UE 2, are installed with the same antennas as UAV relay. Due to the barriers or buildings, there is no direct link between UE 1 and UE 2. For enhancing the reliability of the relay system, SWIPT technique is adopted at UEs, and the PS scheme is considered. To be specific, UEs split the received signal power into two parts, one for EH and the other for ID. We assume that all terminals know the channel state information (CSI) and the amplify-and-forward way is applied at UAV relay.

Fig. 1
figure 1

UAV-enabled two-way relay system with PS-based UEs

With the half-duplex relay, the communication period block T is equally partitioned into two separate phases, i.e., multiple access (MAC) phase and broadcast (BC) phase. In the MAC phase, two UEs transmit their signals to UAV relay. The received signal at UAV relay can be written as

$$\begin{array}{*{20}l} {{\mathbf{y}}_{\mathrm{R}}} = {{\mathbf{H}}_{1}}{{\mathbf{W}}_{1}}{{\mathbf{x}}_{1}} + {{\mathbf{H}}_{2}}{{\mathbf{W}}_{2}}{{\mathbf{x}}_{2}} + {{\mathbf{n}}_{R}}, \end{array} $$
(1)

where \({{\mathbf {H}}_{i}} \in {\mathbb {C}^{N \times N}}\) is the MIMO channel matrix from UE i to UAV relay node, \({{\mathbf {W}}_{i}} \in {\mathbb {C}^{N \times N}}\) represents the beamforming matrix of UE i, which should be under the transmission power constraints tr(WiWiH)=Pi and PiPi, Max, where Pi and Pi, Max denote the realistic transmission power and the limitation of maximum transmission power at UE i, respectively. \({{\mathbf {x}}_{i}} \in {\mathbb {C}^{N \times 1}}\) is the original signal vector with \(\mathbb {E}\left [ {{{\mathbf {x}}_{i}}{{\mathbf {x}}_{i}}^{H}} \right ] = {{\mathbf {I}}_{N}}, {\left (\cdot \right)^{H}}\) is the conjugate transpose operation of one matrix, \({{\mathbf {n}}_{R}} \in {\mathbb {C}^{N \times 1}}\) represents the zero-mean additive white Gaussian noise (AWGN) at UAV relay with covariance matrix \(\sigma _{R}^{2}{{\mathbf {I}}_{N}}\).

In the BC phase, UAV relay forwards its received signals with a beamforming matrix \({{\mathbf {W}}_{R}} \in {\mathbb {C}^{N \times N}}\), which should be under the transmission power constraints tr(WRRWRH)=PR and PRPR, Max, where \({\mathbf {R}} = {{\mathbf {H}}_{1}}{{\mathbf {W}}_{1}}{\mathbf {W}}_{1}^{H}{\mathbf {H}}_{1}^{H} + {{\mathbf {H}}_{2}}{{\mathbf {W}}_{2}}{\mathbf {W}}_{2}^{H}{\mathbf {H}}_{2}^{H} + \sigma _{R}^{2}{{\mathbf {I}}_{N}}, {P_{R}}\) and PR, Max are the realistic transmission power and the limitation of maximum transmission power at UAV relay, respectively. The signal processed by UAV relay is given by

$$\begin{array}{*{20}l} {{\mathbf{x}}_{R}} = {{\mathbf{W}}_{R}}{{\mathbf{y}}_{R}} = {{\mathbf{W}}_{R}}\left({{{\mathbf{H}}_{1}}{{\mathbf{W}}_{1}}{{\mathbf{x}}_{1}} + {{\mathbf{H}}_{2}}{{\mathbf{W}}_{2}}{{\mathbf{x}}_{2}} + {{\mathbf{n}}_{R}}} \right). \end{array} $$
(2)

The received RF signal at UE i is split into two portions with a PS factor, i.e., ρi∈(0,1). ρi portion is used for ID, and 1−ρi portion is for EH. In order to reduce the complexity of MIMO relay system, we assume that PS factors of all antennas at UE i are the same. The ID signal at UE i can be written as

$$\begin{array}{*{20}l} {}{{\mathbf{y}}_{i}} &= \sqrt {{\rho_{i}}} \left({{{\mathbf{G}}_{i}}{{\mathbf{W}}_{R}}{{\mathbf{H}}_{i}}{{\mathbf{W}}_{i}}{{\mathbf{x}}_{i}} + {{\mathbf{G}}_{i}}{{\mathbf{W}}_{R}}{{\mathbf{H}}_{3 - i}}{{\mathbf{W}}_{3 - i}}{{\mathbf{x}}_{3 - i}}}\right.\\ &\qquad\left.{+ {{\mathbf{G}}_{i}}{{\mathbf{W}}_{R}}{{\mathbf{n}}_{R}} + {{\mathbf{n}}_{i}}} \right) + {{\mathbf{n}}_{i,z}}, \end{array} $$
(3)

where ni indicates the zero-mean AWGN vector at UE i with covariance matrix \(\sigma _{i}^{2}{{\mathbf {I}}_{N}}, {{\mathbf {n}}_{i,z}}\) is the additional PS signal processing noise with covariance matrix \(\sigma _{i,z}^{2}{{\mathbf {I}}_{N}}\). 3−i is used to denote the desire signal of paired user. As stated in [23], the self-interference term GiWRHiWixi can be integrally canceled. Therefore, the received signal at UE i can be further obtained as

$$\begin{array}{*{20}l} {}{{\mathbf{y}}_{i}} = \sqrt {{\rho_{i}}} \left({{{\mathbf{G}}_{i}}{{\mathbf{W}}_{R}}{{\mathbf{H}}_{3 - i}}{{\mathbf{W}}_{3 - i}}{{\mathbf{x}}_{3 - i}} + {{\mathbf{G}}_{i}}{{\mathbf{W}}_{R}}{{\mathbf{n}}_{R}} + {{\mathbf{n}}_{i}}} \right) + {{\mathbf{n}}_{i,z}}. \end{array} $$
(4)

According to (4), the received SNR at UE i can be represented by

$$\begin{array}{*{20}l} SN{R_{i}} = \frac{{{{\left\| {{{\mathbf{G}}_{i}}{{\mathbf{W}}_{R}}{{\mathbf{H}}_{3 - i}}{{\mathbf{W}}_{3 - i}}} \right\|}^{2}}}}{{{{\left\| {{{\mathbf{G}}_{i}}{{\mathbf{W}}_{R}}} \right\|}^{2}}\sigma_{R}^{2} + {{\left\| {\sigma_{i}^{2}{{\mathbf{I}}_{N}}} \right\|}^{2}} + {{\left\| {\frac{{\sigma_{i,z}^{2}}}{{{\rho_{i}}}}{{\mathbf{I}}_{N}}} \right\|}^{2}}}}. \end{array} $$
(5)

The harvested energy at UE i can be given by

$$\begin{array}{*{20}l} E_{i,e}^{H} = \frac{T}{2}\eta \left({1 - {\rho_{i}}} \right){\text{tr}}\left({{{\mathbf{G}}_{i}}{{\mathbf{W}}_{R}}{\mathbf{RW}}_{R}^{H}{\mathbf{G}}_{i}^{H} + \sigma_{i}^{2}{{\mathbf{I}}_{N}}} \right), \end{array} $$
(6)

where η∈(0,1) is the energy conversion efficiency of UE i. It is worth pointing out that the self-interference term GiWRHiWixi can also be utilized for EH.

3.2 Problem formulation

In this paper, we design the optimum beamforming matrixes at two UEs and UAV relay and two optimum PS ratios for SWIPT in two-way relay system to maximize the sum retained energy while satisfying sufficient SNR at two UEs and the transmit power constraint at each node. The retained energy of UE i can be expressed as

$$\begin{array}{*{20}l} E_{i}^{R} = E_{i}^{INI}{\text{ + }}E_{i,e}^{H} - \frac{T}{2}{P_{i}}, \end{array} $$
(7)

where \(E_{i}^{INI}\) denotes the initial energy at UE i. The sum retained energy of two UEs can be expressed as

$$\begin{array}{*{20}l} E_{sum}^{R} = \sum\limits_{i = 1}^{2} {E_{i}^{R}}. \end{array} $$
(8)

Based on (5) and (8), the optimization problem of sum retained energy is formulated as

$$\begin{array}{*{20}l} &\underset{\left({{{\mathbf{F}}_1}, {{\mathbf{F}}_2}, {{\mathbf{F}}_R}, {\rho_{1}}, {\rho_{2}}}\right)}{\max} {\text{E}_{sum}^{R}} \end{array} $$
(9a)
$$\begin{array}{*{20}l} &\text{s.t.}\;\;SN{R_{i}} \geqslant {\gamma_{i}},i \in \left\{ {1,2} \right\}, \end{array} $$
(9b)
$$\begin{array}{*{20}l} &\;\;\;\;\;\;\operatorname{Tr} \left({{{\mathbf{W}}_{R}}{\mathbf{RW}}_{R}^{H}} \right) \leqslant {P_{R,\;Max}}, \end{array} $$
(9c)
$$\begin{array}{*{20}l} &\;\;\;\;\;\;\operatorname{Tr} \left({{{\mathbf{W}}_{i}}{{\mathbf{W}}_{i}}^{H}} \right) \leqslant {P_{i,\;Max}},i \in \left\{ {1,2} \right\}, \end{array} $$
(9d)
$$\begin{array}{*{20}l} &\;\;\;\;\;\;\;0 < {\rho_{i}} < 1\;,i \in \left\{ {1,2} \right\}, \end{array} $$
(9e)

where γi indicates the SNR threshold of UE i, constraint (9b) represents the minimum SNR requirement of UE i, (9c) and (9d) denote the maximum transmission power of UAV relay and UE i, respectively, and (9e) is the PS factor constraint of UE i. Since the optimization problem is non-convex, it is intractable to obtain the optimal solution of maximizing the sum retained energy.

4 The proposed alternating optimization scheme

In this section, we propose an AO scheme to convert the complicated optimization problem into three tractable subproblems. Firstly, we optimize the UE beamforming matrixes W1 and W2 with fixed WR,ρ1 and ρ2. Secondly, we try to obtain the optimal relay beamforming matrix WR by assuming the other parameters are given. Finally, the optimal PS factors, i.e., ρ1 and ρ2, are investigated.

4.1 Optimize two UE beamforming matrixes (W1 and W2)

With the above analysis, we first design the optimal solution of W1 and W2 for maximizing the sum retained energy. Substituting \({\mathbf {R}} = {{\mathbf {H}}_{1}}{{\mathbf {W}}_{1}}{\mathbf {W}}_{1}^{H}{{\mathbf {H}}_{1}}^{H} + {{\mathbf {H}}_{2}}{{\mathbf {W}}_{2}}{\mathbf {W}}_{2}^{H}{{\mathbf {H}}_{2}}^{H} + \sigma _{R}^{2}{{\mathbf {I}}_{N}}\) into (9), we have the objective function as

$$ \begin{aligned} E_{sum}^{R} &=\! \frac{T}{2}\;{\text{tr}}\!\left(\!{{\alpha_{1}}{{\mathbf{T}}_{1}}{{\mathbf{W}}_{1}}{\mathbf{W}}_{1}^{H}{\mathbf{T}}_{1}^{H} \,+\, {\alpha_{1}}{{\mathbf{T}}_{3}}{{\mathbf{W}}_{2}}{\mathbf{W}}_{2}^{H}{\mathbf{T}}_{3}^{H} \,-\, {{\mathbf{W}}_{1}}{\mathbf{W}}_{1}^{H}} \right) \hfill \\ & \quad\!\,+\, \frac{T}{2} {\text{tr}}\!\left(\!{{\alpha_{2}}{{\mathbf{T}}_{2}}{{\mathbf{W}}_{1}}{\mathbf{W}}_{1}^{H}{\mathbf{T}}_{2}^{H} \,+\, {\alpha_{2}}{{\mathbf{T}}_{4}}{{\mathbf{W}}_{2}}{\mathbf{W}}_{2}^{H}{\mathbf{T}}_{4}^{H} \,-\, {{\mathbf{W}}_{2}}{\mathbf{W}}_{2}^{H}} \!\right) + J, \end{aligned} $$
(10)

where \({\alpha _{1}} = \eta \left ({1 - {\rho _{1}}} \right), {\alpha _{2}} = \eta \left ({1 - {\rho _{2}}} \right), {{\mathbf {T}}_{1}} = {{\mathbf {G}}_{1}}{{\mathbf {W}}_{R}}{{\mathbf {H}}_{1}}, {{\mathbf {T}}_{2}} = {{\mathbf {G}}_{2}}{{\mathbf {W}}_{R}}{{\mathbf {H}}_{1}}, {{\mathbf {T}}_{3}} = {{\mathbf {G}}_{1}}{{\mathbf {W}}_{R}}{{\mathbf {H}}_{2}}, {{\mathbf {T}}_{4}} = {{\mathbf {G}}_{2}}{{\mathbf {W}}_{R}}{{\mathbf {H}}_{2}}, {{\mathbf {T}}_{5}} = {{\mathbf {G}}_{1}}{{\mathbf {W}}_{R}}{\mathbf {W}}_{R}^{H}{\mathbf {G}}_{1}^{H}, {{\mathbf {T}}_{6}} = {{\mathbf {G}}_{2}}{{\mathbf {W}}_{R}}{\mathbf {W}}_{R}^{H}{\mathbf {G}}_{2}^{H}\), and \(J = \frac {T}{2}{\alpha _{1}}\;{\text {tr}}\left ({\sigma _{R}^{2}{{\mathbf {T}}_{5}} + {{\mathbf {I}}_{N}}\sigma _{1}^{2}} \right) + \frac {T}{2}{\alpha _{2}}\;{\text {tr}}\left ({\sigma _{R}^{2}{{\mathbf {T}}_{6}} + {{\mathbf {I}}_{N}}\sigma _{2}^{2}} \right) + E_{1}^{INI}{\text { + }}E_{2}^{INI}\). Based on the trace property tr(AB)=tr(BA), the total retained energy can be rewritten as

$$\begin{array}{*{20}l} {}E_{sum}^{R} = &\frac{T}{2}\;{\text{tr}}\left({{{\mathbf{\Psi }}_{1}}{{\mathbf{W}}_{1}}{\mathbf{W}}_{1}^{H} + {{\mathbf{\Psi }}_{2}}{{\mathbf{W}}_{1}}{\mathbf{W}}_{1}^{H} - {{\mathbf{W}}_{1}}{\mathbf{W}}_{1}^{H}} \right) \hfill \\ &+ \frac{T}{2}\;{\text{tr}}\left({{{\mathbf{\Psi }}_{3}}{{\mathbf{W}}_{2}}{\mathbf{W}}_{2}^{H} + {{\mathbf{\Psi }}_{4}}{{\mathbf{W}}_{2}}{\mathbf{W}}_{2}^{H} - {{\mathbf{W}}_{2}}{\mathbf{W}}_{2}^{H}} \right) + J \hfill, \end{array} $$
(11)

where \({{\mathbf {\Psi }}_{1}} = {\alpha _{1}}{\mathbf {T}}_{1}^{H}{{\mathbf {T}}_{1}}, {{\mathbf {\Psi }}_{2}} = {\alpha _{2}}{\mathbf {T}}_{2}^{H}{{\mathbf {T}}_{2}}, {{\mathbf {\Psi }}_{3}} = {\alpha _{1}}{\mathbf {T}}_{3}^{H}{{\mathbf {T}}_{3}}\), and \({{\mathbf {\Psi }}_{4}} = {\alpha _{2}}{\mathbf {T}}_{4}^{H}{{\mathbf {T}}_{4}}\). With the trace property in [19], we yield

$$\begin{array}{*{20}l} {\text{tr}}\left({{\mathbf{ABCD}}} \right){\text{ = }}{\left({vec\left({{{\mathbf{D}}^{T}}} \right)} \right)^{T}}\left({{{\mathbf{C}}^{T}} \otimes {\mathbf{A}}} \right)vec\left({\mathbf{B}} \right), \end{array} $$
(12)

where vec(·) is the matrix vectorization operator, ⊗ denotes the Kronecker product. From (12), we can reform the terms \({\text {tr}}\left ({{{\mathbf {\Psi }}_{1}}{{\mathbf {W}}_{1}}{\mathbf {W}}_{1}^{H}} \right)\) and \({\text {tr}}\left ({{{\mathbf {W}}_{1}}{\mathbf {W}}_{1}^{H}} \right)\) into their following forms respectively

$$\begin{array}{*{20}l} {\text{tr}}\left({{{\mathbf{\Psi }}_{1}}{{\mathbf{W}}_{1}}{\mathbf{W}}_{1}^{H}} \right) &= {\text{tr}}\left({{{\mathbf{\Psi }}_{1}}{{\mathbf{W}}_{1}}{{\mathbf{I}}_{N}}{\mathbf{W}}_{1}^{H}} \right) \hfill \\ &= {\mathbf{w}}_{1}^{H}\left({{{\mathbf{I}}_{N}} \otimes {{\mathbf{\Psi }}_{1}}} \right){{\mathbf{w}}_{1}} \hfill \\ &{\text{ = tr}}\left({\left({{{\mathbf{I}}_{N}} \otimes {{\mathbf{\Psi }}_{1}}} \right){{\mathbf{w}}_{1}}{\mathbf{w}}_{1}^{H}} \right) \hfill, \end{array} $$
(13)

and

$$\begin{array}{*{20}l} {\text{tr}}\left({{{\mathbf{W}}_{1}}{\mathbf{W}}_{1}^{H}} \right) = {\mathbf{w}}_{1}^{H}{{\mathbf{w}}_{1}} = {\text{tr}}\left({{{\mathbf{w}}_{1}}{\mathbf{w}}_{1}^{H}} \right), \end{array} $$
(14)

where w1=vec(W1). Similarly, we can convert other terms in (9). We define \({{\mathbf {\Phi }}_{i}} \triangleq {{\mathbf {w}}_{i}}{\mathbf {w}}_{i}^{H}\), then the sum retained energy becomes

$$\begin{array}{*{20}l} {}E_{sum}^{R}= &\frac{T}{2}\;{\text{tr}}\left({\left({\left({{{\mathbf{I}}_{N}} \otimes {{\mathbf{\Psi }}_{1}}} \right) + \left({{{\mathbf{I}}_{N}} \otimes {{\mathbf{\Psi }}_{2}}} \right) - {{\mathbf{I}}_{E}}} \right){{\mathbf{\Phi }}_{1}}} \right) \hfill \\ & + \frac{T}{2}\;{\text{tr}}\left({\left({\left({{{\mathbf{I}}_{N}} \otimes {{\mathbf{\Psi }}_{3}}} \right) + \left({{{\mathbf{I}}_{N}} \otimes {{\mathbf{\Psi }}_{4}}} \right) - {{\mathbf{I}}_{E}}} \right){{\mathbf{\Phi }}_{2}}} \right) + J \hfill, \end{array} $$
(15)

where \({{\mathbf {I}}_{E}} \in {\mathbb {C}^{{N^{2}} \times {N^{2}}}}\) is a unit matrix. With (5), we can transform the UEs SNR constraints into

$$\begin{array}{*{20}l} {\text{tr}}\left({\left({{{\mathbf{I}}_{N}} \otimes {\mathbf{T}}_{3}^{H}{{\mathbf{T}}_{3}}} \right){{\mathbf{\Phi }}_{2}}} \right) \geqslant {\gamma_{1}}{J_{1}}, \end{array} $$
(16)
$$\begin{array}{*{20}l} {\text{tr}}\left({\left({{{\mathbf{I}}_{N}} \otimes {\mathbf{T}}_{2}^{H}{{\mathbf{T}}_{2}}} \right){{\mathbf{\Phi }}_{1}}} \right) \geqslant {\gamma_{2}}{J_{2}}, \end{array} $$
(17)

where \({J_{1}} = \left \| {{{\mathbf {G}}_{1}}{{\mathbf {F}}_{R}}} \right \|_{}^{2}\sigma _{R}^{2} + \sigma _{1}^{2} + \frac {{\sigma _{1,z}^{2}}}{{{\rho _{1}}}}\) and \({J_{2}} = \left \| {{{\mathbf {G}}_{2}}{{\mathbf {F}}_{R}}} \right \|_{}^{2}\sigma _{R}^{2} + \sigma _{2}^{2} + \frac {{\sigma _{2,z}^{2}}}{{{\rho _{2}}}}\). Similarly, the UAV relay transmission power constraint will become

$$\begin{array}{*{20}l} {}{\text{tr}}\left({\left({{{\mathbf{I}}_{N}} \otimes {{\mathbf{T}}_{5}}} \right){{\mathbf{\Phi }}_{1}} + \left({{{\mathbf{I}}_{N}} \otimes {{\mathbf{T}}_{6}}} \right){{\mathbf{\Phi }}_{2}} + \sigma_{R}^{2}{{\mathbf{W}}_{R}}{\mathbf{W}}_{R}^{H}} \right) \leqslant {P_{R,\;Max}}. \end{array} $$
(18)

Plugging (16), (17), and (18) back into (9), the original optimization problem can be reformulated as

$$\begin{array}{*{20}l} &{} \underset{\left({{{\mathbf{\Phi }}_{1}},\;{{\mathbf{\Phi }}_{2}} \succcurlyeq {\mathbf{0}}} \right)}{\max }\;\;\;\frac{T}{2}\;{\text{tr}}\left({{{\mathbf{Q}}_{1}}{{\mathbf{\Phi }}_{1}}} \right) + \frac{T}{2}\;{\text{tr}}\left({{{\mathbf{Q}}_{2}}{{\mathbf{\Phi }}_{2}}} \right) + {J} \end{array} $$
(19a)
$$\begin{array}{*{20}l} &{} s.t.\;\;{\text{tr}}\left({\left({{{\mathbf{I}}_{N}} \otimes {\mathbf{T}}_{3}^{H}{{\mathbf{T}}_{3}}} \right){{\mathbf{\Phi }}_{2}}} \right) \geqslant {\gamma_{1}}{J_{1}}, \end{array} $$
(19b)
$$\begin{array}{*{20}l} &{} \;\;\;\;\;\;\;{\text{tr}}\left({\left({{{\mathbf{I}}_{N}} \otimes {\mathbf{T}}_{2}^{H}{{\mathbf{T}}_{2}}} \right){{\mathbf{\Phi }}_{1}}} \right) \geqslant {\gamma_{2}}{J_{2}} \hfill, \end{array} $$
(19c)
$$\begin{array}{*{20}l} & {} \;\;\;\;\;\;\;{\text{tr}}\left({{{\mathbf{\Phi }}_{i}}} \right) \leqslant {P_{i,\;Max}},\;i \in \left\{ {1,2} \right\}, \end{array} $$
(19d)
$$\begin{array}{*{20}l} & {} \;\;\;\;\;\;\;{\text{tr}}\left({\left({{{\mathbf{I}}_{N}} \otimes {{\mathbf{T}}_{5}}} \right){{\mathbf{\Phi }}_{1}} + \left({{{\mathbf{I}}_{N}} \otimes {{\mathbf{T}}_{6}}} \right){{\mathbf{\Phi }}_{2}} + \sigma_{R}^{2}{{\mathbf{W}}_{R}}{\mathbf{W}}_{R}^{H}} \right)\\ {}&\qquad\leqslant {P_{R,\;Max}}, \end{array} $$
(19e)

where Q1 = (INΨ1)+(INΨ2)−IE and Q2 = (INΨ3)+(INΨ4)−IE. According to [27], the optimization problem (19) can be converted into a semi-definite programming (SDP) problem without a rank-one constraint. Due to the linearity of the objective function and the relaxed constraints, problem (18) which is convex can be solved by the classical optimization tools, e.g., CVX. Therefore, with the optimal solutions \({\mathbf {\Phi }}_{1}^{\text {*}}\) and \({\mathbf {\Phi }}_{2}^{\text {*}}\), we can obtain the optimal beamforming matrixes of UEs, i.e., \({\mathbf {W}}_{1}^{\text {*}}\) and \({\mathbf {W}}_{2}^{\text {*}}\).

4.2 Optimize UAV relay beamforming (WR)

Similarly, with given W1,W2 and two UE PS factors, i.e., ρ1,ρ2, the total retained energy of two UEs is equivalent to

$$\begin{array}{*{20}l} E_{sum}^{R} = \frac{T}{2}{\text{tr}}\left({{{{\tilde{\mathbf{{\Psi}}}}}_{1}}{{\mathbf{W}}_{R}}{\mathbf{RW}}_{R}^{H} + {{{\tilde{\mathbf{{\Psi}}}}}_{2}}{{\mathbf{W}}_{R}}{\mathbf{RW}}_{R}^{H}} \right) + \tilde{J}, \end{array} $$
(20)

where \({{\tilde {\mathbf {\Psi }}}_{1}} = {\alpha _{1}}{\mathbf {G}}_{1}^{H}{{\mathbf {G}}_{1}}, {{\tilde {\mathbf {\Psi }}}_{2}} = {\alpha _{2}}{\mathbf {G}}_{2}^{H}{{\mathbf {G}}_{2}}\), and \(\tilde {J} = \frac {T}{2}{\text {tr}}\)\(\left ({{\alpha _{1}}\sigma _{1}^{2}{{\mathbf {I}}_{N}} + {\alpha _{2}}\sigma _{2}^{2}{{\mathbf {I}}_{N}}} \right) + E_{1}^{INI}{\text { + }}E_{2}^{INI} - \frac {T}{2}\left ({{P_{1}} + {P_{2}}} \right)\). Based on (12), we have

$$\begin{array}{*{20}l} {}E_{sum}^{R} = \frac{T}{2}{\text{tr}}\left({\left({\left({{{\mathbf{R}}^{T}} \otimes {{{\tilde{\mathbf{{\Psi}}}}}_{1}}} \right) + \left({{{\mathbf{R}}^{T}} \otimes {{{\tilde{\mathbf{{\Psi}}}}}_{2}}} \right)} \right){{\mathbf{w}}_{R}}{\mathbf{w}}_{R}^{H}} \right) + \tilde{J}, \end{array} $$
(21)

where wR=vec(WR). Further, the constraints in (9) can be rewritten as

$$\begin{array}{*{20}l} \text{tr}\left({{{\mathbf{K}}_{1}}{{\mathbf{w}}_{R}}{\mathbf{w}}_{R}^{H}} \right) \geqslant {\gamma_{1}}\left({\sigma_{1}^{2} + \frac{{\sigma_{1,z}^{2}}}{{{\rho_{1}}}}} \right), \end{array} $$
(22)
$$\begin{array}{*{20}l} \text{tr}\left({{{\mathbf{K}}_{2}}{{\mathbf{w}}_{R}}{\mathbf{w}}_{R}^{H}} \right) \geqslant {\gamma_{2}}\left({\sigma_{2}^{2} + \frac{{\sigma_{2,z}^{2}}}{{{\rho_{2}}}}} \right), \end{array} $$
(23)

and

$$\begin{array}{*{20}l} \text{tr}\left({\left({{{\mathbf{R}}^{T}} \otimes {{\mathbf{I}}_{N}}} \right){{\mathbf{w}}_{R}}{\mathbf{w}}_{R}^{H}} \right) \leqslant {P_{R,\;Max}}, \end{array} $$
(24)

where

$$\begin{array}{*{20}l} {}{{\mathbf{K}}_{1}}{\text{ = }}\left({{{\left({{{\mathbf{H}}_{2}}{{\mathbf{F}}_{2}}{\mathbf{F}}_{2}^{H}{\mathbf{H}}_{2}^{H}} \right)}^{T}} \otimes \left({{\mathbf{G}}_{1}^{H}{{\mathbf{G}}_{1}}} \right)} \right){\text{ - }}{\gamma_{1}}\sigma_{R}^{2}\left({{{\mathbf{I}}_{N}} \otimes \left({{\mathbf{G}}_{1}^{H}{{\mathbf{G}}_{1}}} \right)} \right), \end{array} $$
(25)
$$\begin{array}{*{20}l} {}{{\mathbf{K}}_{2}}{\text{ = }}\left({{{\left({{{\mathbf{H}}_{1}}{{\mathbf{F}}_{1}}{\mathbf{F}}_{1}^{H}{\mathbf{H}}_{1}^{H}} \right)}^{T}} \otimes \left({{\mathbf{G}}_{2}^{H}{{\mathbf{G}}_{2}}} \right)} \right){\text{ - }}{\gamma_{2}}\sigma_{R}^{2}\left({{{\mathbf{I}}_{N}} \otimes \left({{\mathbf{G}}_{2}^{H}{{\mathbf{G}}_{2}}} \right)} \right). \end{array} $$
(26)

Here, we define \({{\mathbf {\Phi }}_{R}} \triangleq {{\mathbf {w}}_{R}}{\mathbf {w}}_{R}^{H}\). Substituting (22), (23), and (24) back into (9), the original optimization problem can be reformulated as

$$\begin{array}{*{20}l} &{} \underset{\left({{{\mathbf{\Phi }}_{R}} \succcurlyeq {\mathbf{0}}} \right)} {\max}\;\frac{T}{2}{\text{tr}}\left({\left({\left({{{\mathbf{R}}^{T}} \otimes {{{\tilde{\mathbf{{\Psi}}}}}_{1}}} \right) + \left({{{\mathbf{R}}^{T}} \otimes {{{\tilde{\mathbf{{\Psi}}}}}_{2}}} \right)} \right){{\mathbf{\Phi }}_{R}}} \right) + {\tilde{J}} \end{array} $$
(27a)
$$\begin{array}{*{20}l} &{} s.t.\;\;\text{tr}\left({{{\mathbf{K}}_{1}}{{\mathbf{\Phi }}_{R}}} \right) \geqslant {\gamma_{1}}\left({\sigma_{1}^{2} + \frac{{\sigma_{1,z}^{2}}}{{{\rho_{1}}}}} \right) \hfill, \end{array} $$
(27b)
$$\begin{array}{*{20}l} & {}\;\;\;\;\;\;\; \text{tr}\left({{{\mathbf{K}}_{2}}{{\mathbf{\Phi }}_{R}}} \right) \geqslant {\gamma_{2}}\left({\sigma_{2}^{2} + \frac{{\sigma_{2,z}^{2}}}{{{\rho_{2}}}}} \right) \hfill, \end{array} $$
(27c)
$$\begin{array}{*{20}l} &{} \;\;\;\;\;\;\; \text{tr}\left({\left({{{\mathbf{R}}^{T}} \otimes {{\mathbf{I}}_{N}}} \right){{\mathbf{\Phi }}_{R}}} \right) \leqslant {P_{R,\;Max}}. \end{array} $$
(27d)

Similar to the first subproblem, we can utilize SDP scheme to solve the problem (27). With the optimal solutions \({\mathbf {\Phi }}_{R}^{\text {*}}\), we can get the optimal beamforming matrixes of UAV relay, i.e., \({\mathbf {W}}_{R}^{\text {*}}\).

4.3 Optimize two PS factors (ρ1,ρ1)

For given W1,W2, and WR, the objective function of maximizing the sum retained energy with ρ1 and ρ2 can be written as

$$\begin{array}{*{20}l} E_{sum}^{R} = {\beta_{1}}\left({1 - {\rho_{1}}} \right){\text{ + }}{\beta_{2}}\left({1 - {\rho_{2}}} \right){\text{ + }}{\hat{J}}, \end{array} $$
(28)

where \({\beta _{1}} = \frac {T}{2}\eta {\text {tr}}\left ({{{\mathbf {G}}_{1}}{{\mathbf {W}}_{R}}{\mathbf {RW}}_{R}^{H}{\mathbf {G}}_{1}^{H} + \sigma _{1}^{2}{{\mathbf {I}}_{N}}} \right), {\beta _{2}} = \frac {T}{2}\eta {\text {tr}} \left ({{{\mathbf {G}}_{2}}{{\mathbf {W}}_{R}}{\mathbf {RW}}_{R}^{H}{\mathbf {G}}_{2}^{H} \,+\, \sigma _{2}^{2}{{\mathbf {I}}_{N}}} \!\right)\), and \({\hat {J}} \,=\, E_{1}^{INI}{\text { \,+\, }}E_{2}^{INI} - \frac {T}{2}\left ({{P_{1}} + {P_{2}}} \right)\). According to the SNR constraints in (9), we obtain

$$\begin{array}{*{20}l} {\rho_{1}} \geqslant \frac{{\sigma_{1,z}^{2}{\gamma_{1}}}}{{{\mu_{1}}}}, \end{array} $$
(29)
$$\begin{array}{*{20}l} {\rho_{2}} \geqslant \frac{{\sigma_{2,z}^{2}{\gamma_{2}}}}{{{\mu_{2}}}}, \end{array} $$
(30)

where \({\mu _{1}} = {\left \| {{{\mathbf {G}}_{1}}{{\mathbf {W}}_{R}}{{\mathbf {H}}_{2}}{{\mathbf {W}}_{2}}} \right \|^{2}} - {\gamma _{1}}\sigma _{R}^{2}{\left \| {{{\mathbf {G}}_{1}}{{\mathbf {W}}_{R}}} \right \|^{2}} - {\gamma _{1}}\sigma _{1}^{2}\) and \({\mu _{2}} = {\left \| {{{\mathbf {G}}_{2}}{{\mathbf {W}}_{R}}{{\mathbf {H}}_{1}}{{\mathbf {F}}_{1}}} \right \|^{2}} - {\gamma _{2}}\sigma _{R}^{2}{\left \| {{{\mathbf {G}}_{2}}{{\mathbf {W}}_{R}}} \right \|^{2}} - {\gamma _{2}}\sigma _{2}^{2}\). Therefore, the original optimization problem can be changed as

$$\begin{array}{*{20}l} &\underset{\left({{\rho_{1}},{\rho_{2}}} \right)}{\max} \;{\beta_{1}}\left({1 - {\rho_{1}}} \right){\text{ + }}{\beta_{2}}\left({1 - {\rho_{2}}} \right){\text{ + }}{\hat {J}} \hfill \end{array} $$
(31a)
$$\begin{array}{*{20}l} & s.t.\;\;\;{\rho_{1}} \geqslant \frac{{\sigma_{1,z}^{2}{\gamma_{1}}}}{{{\mu_{1}}}} \hfill, \end{array} $$
(31b)
$$\begin{array}{*{20}l} &\;\;\;\;\;\;\;\;{\rho_{2}} \geqslant \frac{{\sigma_{2,z}^{2}{\gamma_{2}}}}{{{\mu_{2}}}} \hfill, \end{array} $$
(31c)
$$\begin{array}{*{20}l} & \;\;\;\;\;\;\;\;0 < {\rho_{1}} < 1 \hfill, \end{array} $$
(31d)
$$\begin{array}{*{20}l} & \;\;\;\;\;\;\;\;0 < {\rho_{2}} < 1 \hfill. \end{array} $$
(31e)

With the known feasible region of the problem (31) and the relationship of the variates, the optimal closed-form expression of PS factors can be expressed as

$$\begin{array}{*{20}l} {\rho_{1,\;opt}} = \frac{{\sigma_{1,z}^{2}{\gamma_{1}}}}{{{\mu_{1}}}}, \end{array} $$
(32)
$$\begin{array}{*{20}l} {\rho_{2,\;opt}} = \frac{{\sigma_{2,z}^{2}{\gamma_{2}}}}{{{\mu_{2}}}}. \end{array} $$
(33)

In summary, the proposed AO scheme for maximizing the sum remained energy can be summarized as Algorithm 1.

Complexity Analysis: It is seen that the optimization problems (19) and (27) are two classical SDP problems which can be solved by utilizing the interior-point method (IPM). According to [18], the complexity of the proposed AO scheme for maximizing the sum retained energy is about \(\mathcal {O}\left ({{M_{Iter}} \cdot 3{N^{7}}} \right)\), where MIter indicates the number of iterations.

5 Low-complexity beamforming design based on GSVD technique

The proposed AO scheme of UAV-enable two-way relay system in Section 4 gets great rate performance, whereas it has high computational complexity. To provide an outstanding balance between complexity and performance, a low-complexity scheme based on GSVD technique is presented in this section. Compared with the above AO scheme, the low-complexity GSVD scheme transforms the original optimization problem with matrix variables into the power allocation optimization problem with scalar variables. However, two schemes have the similar iteration process.

5.1 Beamforming design based on GSVD

By adopting GSVD technique onto the related MAC channel \(\left [ {{\mathbf {H}}_{1}^{H},\;{\mathbf {H}}_{2}^{H}} \right ]\), we can decompose H1 and H2 as

$$\begin{array}{*{20}l} {{\mathbf{H}}_{1}} = {{\mathbf{B}}_{h}}{{\mathbf{\Sigma }}_{1}}{\mathbf{A}}_{{h_{1}}}^{H},\;\;\;\;{{\mathbf{H}}_{2}} = {{\mathbf{B}}_{h}}{{\mathbf{\Sigma }}_{2}}{\mathbf{A}}_{{h_{2}}}^{H}, \end{array} $$
(34)

where Bh denotes a full-rank complex matrix, \({{\mathbf {\Sigma }}_{1}} \in {\mathbb {C}^{N \times N}}\) and \({{\mathbf {\Sigma }}_{2}} \in {\mathbb {C}^{N \times N}}\) are two diagonal matrixes, \({{\mathbf {A}}_{{h_{1}}}} \in {\mathbb {C}^{N \times N}}\) and \({{\mathbf {A}}_{{h_{2}}}} \in {\mathbb {C}^{N \times N}}\) represent two unitary matrixes. After processing the MAC channel, we apply SVD technique to decompose the BC channel \({{\mathbf {G}}_{BC}} = {\left [ {{\mathbf {G}}_{1}^{T},\;{\mathbf {G}}_{2}^{T}} \right ]^{T}}\) as

$$\begin{array}{*{20}l} {{\mathbf{G}}_{BC}} = {{\mathbf{B}}_{g}}{{\mathbf{\Sigma }}_{g}}{\mathbf{A}}_{g}^{H}, \end{array} $$
(35)

where \({{\mathbf {B}}_{g}} \in {\mathbb {C}^{2N \times 2N}}\) and \({{\mathbf {A}}_{g}} \in {\mathbb {C}^{N \times N}}\) denotes two unitary matrixes, and \({{\mathbf {\Sigma }}_{g}} \in {\mathbb {C}^{N \times N}}\) is a diagonal matrix with no-negative real values arranged into decreasing order. Based on (35), G1 and G2 are decomposed as [17, 28]

$$\begin{array}{*{20}l} {{\mathbf{G}}_{1}} = {{\mathbf{B}}_{{g_{1}}}}{{\mathbf{\Sigma }}_{g}}{\mathbf{A}}_{g}^{H},\;\;\;{{\mathbf{G}}_{2}} = {{\mathbf{B}}_{{g_{2}}}}{{\mathbf{\Sigma }}_{g}}{\mathbf{A}}_{g}^{H}, \end{array} $$
(36)

where \({{\mathbf {B}}_{{g_{1}}}} = {{\mathbf {B}}_{g}}\left ({1:N,1:N} \right)\) and \({{\mathbf {B}}_{{g_{2}}}} = {{\mathbf {B}}_{g}}\left ({N + 1:2N,1:N} \right)\). It is worth noting that \({{\mathbf {B}}_{{g_{i}}}}\) is not a unitary matrix. According to the above decompositions of MAC and BC channel matrixes, the beamforming structures of UE i and UAV relay can be proposed as

$$\begin{array}{*{20}l} {{\mathbf{W}}_{i}} = {{\mathbf{A}}_{{h_{i}}}}{{\mathbf{\Lambda }}_{i}},\;\;\;{{\mathbf{W}}_{R}} = {{\mathbf{A}}_{g}}{{\mathbf{\Lambda }}_{R}}{\mathbf{B}}_{h}^{- 1}, \end{array} $$
(37)

where \({{\mathbf {\Lambda }}_{i}} \in {\mathbb {C}^{N \times N}}\) and \({{\mathbf {\Lambda }}_{R}} \in {\mathbb {C}^{N \times N}}\) are two diagonal matrixes standing for the power allocation of UE i and UAV relay, respectively. Plugging (34), (36), and (37) into (4), the receive signal at UE i can be written as

$$\begin{array}{*{20}l} {}{{\mathbf{y}}_{i}} &= \sqrt {{\rho_{i}}} \left({{{\mathbf{B}}_{{g_{i}}}}{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}{{\mathbf{\Sigma }}_{3 - i}}{{\mathbf{\Lambda }}_{3 - i}}{{\mathbf{x}}_{3 - i}} \,+\, {{\mathbf{B}}_{{g_{i}}}}{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}{\mathbf{B}}_{h}^{- 1}{{\mathbf{n}}_{R}} + {{\mathbf{n}}_{i}}} \right)\\{} &\quad+ {{\mathbf{n}}_{i,z}}. \end{array} $$
(38)

After using the zero-forcing operation for the received signal at UE i, i.e., \({\mathbf {B}}_{{g_{i}}}^{- 1}{{\mathbf {y}}_{i}}\), we have

$$\begin{array}{*{20}l} {}{{\tilde{\mathbf{y}}}_{i}} = \sqrt {{\rho_{i}}} \left({{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}{{\mathbf{\Sigma }}_{3 - i}}{{\mathbf{\Lambda }}_{3 - i}}{{\mathbf{x}}_{3 - i}} + {{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}{{{\tilde{\mathbf{n}}}}_{R}} + {{{\tilde{\mathbf{n}}}}_{i}}} \right) + {{\tilde{\mathbf{n}}}_{i,z}}, \end{array} $$
(39)

where \({{\tilde {\mathbf {y}}}_{i}} = {\mathbf {B}}_{{g_{i}}}^{- 1}{{\mathbf {y}}_{i}}, {{\tilde {\mathbf {n}}}_{R}} = {\mathbf {B}}_{h}^{- 1}{{\mathbf {n}}_{R}}, {{\tilde {\mathbf {n}}}_{i}} = {\mathbf {B}}_{{g_{i}}}^{- 1}{{\mathbf {n}}_{i}}\), and \({{\tilde {\mathbf {n}}}_{i,z}} = {\mathbf {B}}_{{g_{i}}}^{- 1}{{\mathbf {n}}_{i,z}}\). Substituting (37) back into (7), we can transform the retained energy \(E_{i}^{R}\) into its equivalent form as

$$\begin{array}{*{20}l} {}E_{i}^{R} &\,=\, \frac{T}{2}\eta \left({1 \,-\, {\rho_{i}}} \right){\text{tr}}\!\left(\!{{{\mathbf{\Psi }}_{{g_{i}}}}{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}\left({{\mathbf{\Pi }} \,+\, \sigma_{R}^{2}{{\mathbf{\Psi }}_{h}}} \right)\!{\mathbf{\Lambda }}_{R}^{H}{\mathbf{\Sigma }}_{g}^{H} + \sigma_{i}^{2}{{\mathbf{I}}_{N}}} \right) \hfill \\ &\quad{\text{ + }}E_{i}^{INI} - \frac{T}{2}{\text{tr}}\left({{{\mathbf{\Lambda }}_{i}}{\mathbf{\Lambda }}_{i}^{H}} \right) \hfill, \end{array} $$
(40)

where \({{\mathbf {\Psi }}_{{g_{i}}}} = {\mathbf {B}}_{{g_{i}}}^{H}{{\mathbf {B}}_{{g_{i}}}}, {{\mathbf {\Psi }}_{h}} = {\left ({{\mathbf {B}}_{h}^{H}{{\mathbf {B}}_{h}}} \right)^{- 1}}\), and \({\mathbf {\Pi }} = {{\mathbf {\Sigma }}_{1}}{{\mathbf {\Lambda }}_{1}}{\mathbf {\Lambda }}_{1}^{H}{\mathbf {\Sigma }}_{1}^{H} + {{\mathbf {\Sigma }}_{2}}{{\mathbf {\Lambda }}_{2}}{\mathbf {\Lambda }}_{2}^{H}{\mathbf {\Sigma }}_{2}^{H}\). With the transformation of the optimization constraints, the problem (9) can be reformulated as

$$\begin{array}{*{20}l} &{} \underset{\left({{{\mathbf{\Lambda }}_{1}},\;{{\mathbf{\Lambda }}_{2}},\;{{\mathbf{\Lambda }}_{R}},\;{\rho_{1}},\;{\rho_{2}}} \right)}{\max} \;\;\;E_{1}^{R} + E_{2}^{R} \hfill \end{array} $$
(41a)
$$\begin{array}{*{20}l} &{}s.t.\;\;\frac{{{\text{tr}}\left({{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}{{\mathbf{\Sigma }}_{3 - i}}{{\mathbf{\Lambda }}_{3 - i}}{\mathbf{\Lambda }}_{3 - i}^{H}{\mathbf{\Sigma }}_{3 - i}^{H}{\mathbf{\Lambda }}_{R}^{H}{\mathbf{\Sigma }}_{g}^{H}} \right)}}{{{\text{tr}}\left({{\mathbf{\Lambda }}_{R}^{H}{\mathbf{\Sigma }}_{g}^{H}{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}\;{{\mathbf{\Psi }}_{h}}} \right)\sigma_{R}^{2}{\text{ + tr}}\left({{\mathbf{\Psi }}_{{g_{i}}}^{- 1}} \right)\sigma_{i}^{2}{\text{ + tr}}\left({{\mathbf{\Psi }}_{{g_{i}}}^{- 1}} \right)\frac{{\sigma_{i,z}^{2}}}{{{\rho_{i}}}}}} \geqslant {\gamma_{i}},i \in \left\{ {1,2} \right\}, \hfill \end{array} $$
(41b)
$$\begin{array}{*{20}l} &{} \;\;\;\;\;\;{\text{tr}}\left({{{\mathbf{\Lambda }}_{R}}\left({{\mathbf{\Pi }} + {{\mathbf{\Psi }}_{h}}\sigma_{R}^{2}} \right){\mathbf{\Lambda }}_{R}^{H}} \right) \leqslant {P_{R,\;Max}}, \hfill \end{array} $$
(41c)
$$\begin{array}{*{20}l} &{} \;\;\;\;\;\;{\text{tr}}\left({{{\mathbf{\Lambda }}_{i}}{\mathbf{\Lambda }}_{i}^{H}} \right) \leqslant {P_{i,\;Max}},i \in \left\{ {1,2} \right\}, \hfill \end{array} $$
(41d)
$$\begin{array}{*{20}l} &{} \;\;\;\;\;\;0 < {\rho_{i}} < 1\;,i \in \left\{ {1,2} \right\}. \hfill \end{array} $$
(41e)

Obviously, the trace matrixes of the retained energy \(E_{i}^{R}\) and two constraints in (41) are not diagonal. According to Property 1 in [17], we can achieve a lower bound on the retained energy to simplify the optimization problem (41). A lower bound of \(E_{i}^{R}\) is expressed as

$$\begin{array}{*{20}l} {} E_{i}^{R} \!\geqslant \!\overset{{\smash{\scriptscriptstyle\frown}}}{E}_{i}^{R} &\,=\, \frac{T}{2}\eta \!\left({1 \,-\, {\rho_{i}}} \right){\text{tr}}\!\left(\! {{{\mathbf{\Lambda }}_{{g_{i}}}}{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}\!\left(\! {{\mathbf{\Pi }} \!+ \!\sigma_{R}^{2}{{\mathbf{\Lambda }}_{h}}} \!\right)\!{\mathbf{\Lambda }}_{R}^{H}{\mathbf{\Sigma }}_{g}^{H} \,+\, \sigma_{i}^{2}{{\mathbf{I}}_{N}}} \!\right) \\ &\quad {\text{ + }}E_{i}^{INI} - \frac{T}{2}{\text{tr}}\left({{{\mathbf{\Lambda }}_{i}}{\mathbf{\Lambda }}_{i}^{H}} \right), \end{array} $$
(42)

where \({{\mathbf {\Lambda }}_{{g_{i}}}}\) and Λh are two diagonal matrixes that consist of N diagonal elements of \({{\mathbf {\Psi }}_{{g_{i}}}}\) and Ψh. For arbitrary positive semi-definite matrixes A and X, there are \({\text {tr}}\left ({{\mathbf {AX}}} \right) \geqslant {\text {tr}}\left ({{\mathbf {A}}{{\mathbf {\Lambda }}_{X}}} \right)\) in [17] and \({\text {tr}}\left ({{{\mathbf {X}}^{- 1}}} \right) \geqslant \sum \limits _{i} {{{\left [ {{\mathbf {X}}\left ({i,i} \right)} \right ]}^{- 1}}} \) in [26], where ΛX is a diagonal matrix that comprises the diagonal elements of X. Hence, we can obtain the upper bound of SNR constraint (41a) and the lower bound of UAV relay transmission power constraint (41b) respectively as

$$\begin{array}{*{20}l} &{}\frac{{{\text{tr}}\left({{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}{{\mathbf{\Sigma }}_{3 - i}}{{\mathbf{\Lambda }}_{3 - i}}{\mathbf{\Lambda }}_{3 - i}^{H}{\mathbf{\Sigma }}_{3 - i}^{H}{\mathbf{\Lambda }}_{R}^{H}{\mathbf{\Sigma }}_{g}^{H}} \right)}}{{{\text{tr}}\left({{\mathbf{\Lambda }}_{R}^{H}{\mathbf{\Sigma }}_{g}^{H}{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}\;{{\mathbf{\Psi }}_{h}}} \right)\sigma_{R}^{2}{\text{ + tr}}\left({{\mathbf{\Psi }}_{{g_{i}}}^{- 1}} \right)\sigma_{i}^{2}{\text{ + tr}}\left({{\mathbf{\Psi }}_{{g_{i}}}^{- 1}} \right)\frac{{\sigma_{i,z}^{2}}}{{{\rho_{i}}}}}} \hfill \\ &{}\;\; \leqslant \frac{{{\text{tr}}\left({{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}{{\mathbf{\Sigma }}_{3 - i}}{{\mathbf{\Lambda }}_{3 - i}}{\mathbf{\Lambda }}_{3 - i}^{H}{\mathbf{\Sigma }}_{3 - i}^{H}{\mathbf{\Lambda }}_{R}^{H}{\mathbf{\Sigma }}_{g}^{H}} \right)}}{{{\text{tr}}\left({{\mathbf{\Lambda }}_{R}^{H}{\mathbf{\Sigma }}_{g}^{H}{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}\;{{\mathbf{\Lambda }}_{h}}} \right)\sigma_{R}^{2}{\text{ + tr}}\left({{\mathbf{\Lambda }}_{{g_{i}}}^{- 1}} \right)\sigma_{i}^{2}{\text{ + tr}}\left({{\mathbf{\Lambda }}_{{g_{i}}}^{- 1}} \right)\frac{{\sigma_{i,z}^{2}}}{{{\rho_{i}}}}}} \hfill, \end{array} $$
(43)
$$\begin{array}{*{20}l} {}{\text{tr}}\left({{{\mathbf{\Lambda }}_{R}}\left({{\mathbf{\Pi }} + {{\mathbf{\Psi }}_{h}}\sigma_{R}^{2}} \right){\mathbf{\Lambda }}_{R}^{H}} \right) \geqslant {\text{tr}}\left({{{\mathbf{\Lambda }}_{R}}\left({{\mathbf{\Pi }} + {{\mathbf{\Lambda }}_{h}}\sigma_{R}^{2}} \right){\mathbf{\Lambda }}_{R}^{H}} \right). \end{array} $$
(44)

With (42), (43), and (44), we can reformulate the problem (41) as

$$\begin{array}{*{20}l} &{} \underset{\left({{{\mathbf{\Lambda }}_{1}},\;{{\mathbf{\Lambda }}_{2}},\;{{\mathbf{\Lambda }}_{R}},\;{\rho_{1}},\;{\rho_{2}}} \right)}{\max} \;\;\;\overset{{\smash{\scriptscriptstyle\frown}}}{E}_{1}^{R} + \overset{{\smash{\scriptscriptstyle\frown}}}{E}_{2}^{R} \end{array} $$
(45a)
$$\begin{array}{*{20}l} &{} s.t.\;\frac{{{\text{tr}}\left({{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}{{\mathbf{\Sigma }}_{3 - i}}{{\mathbf{\Lambda }}_{3 - i}}{\mathbf{\Lambda }}_{3 - i}^{H}{\mathbf{\Sigma }}_{3 - i}^{H}{\mathbf{\Lambda }}_{R}^{H}{\mathbf{\Sigma }}_{g}^{H}} \right)}}{{{\text{tr}}\left({{\mathbf{\Lambda }}_{R}^{H}{\mathbf{\Sigma }}_{g}^{H}{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}\;{{\mathbf{\Lambda }}_{h}}} \right)\sigma_{R}^{2}{\text{ + tr}}\left({{\mathbf{\Lambda }}_{{g_{i}}}^{- 1}} \right)\sigma_{i}^{2}{\text{ + tr}}\left({{\mathbf{\Lambda }}_{{g_{i}}}^{- 1}} \right)\frac{{\sigma_{i,z}^{2}}}{{{\rho_{i}}}}}}\\ &{} \geqslant {\gamma_{i}},i \in \left\{ {1,2} \right\}, \hfill \end{array} $$
(45b)
$$\begin{array}{*{20}l} &{} \;\;\;\;\;\;{\text{tr}}\left({{{\mathbf{\Lambda }}_{R}}\left({{\mathbf{\Pi }} + {{\mathbf{\Lambda }}_{h}}\sigma_{R}^{2}} \right){\mathbf{\Lambda }}_{R}^{H}} \right) \leqslant {P_{R,\;Max}}, \hfill \end{array} $$
(45c)
$$\begin{array}{*{20}l} &{} \;\;\;\;\;\;{\text{tr}}\left({{{\mathbf{\Lambda }}_{i}}{\mathbf{\Lambda }}_{i}^{H}} \right) \leqslant {P_{i,\;Max}},i \in \left\{ {1,2} \right\}, \hfill \end{array} $$
(45d)
$$\begin{array}{*{20}l} &{} \;\;\;\;\;\;0 < {\rho_{i}} < 1\;,i \in \left\{ {1,2} \right\} \hfill. \end{array} $$
(45e)

5.2 Optimize UE power allocation

In the following, we adopt the AO-based optimization algorithm to solve the intractable and non-convex problem (45). For given ΛR,ρ1, and ρ2, the objective function of the sum retained energy can be rewritten as

$$\begin{array}{*{20}l} {}\overset{{\smash{\scriptscriptstyle\frown}}}{E}_{1}^{R} + \overset{{\smash{\scriptscriptstyle\frown}}}{E}_{2}^{R} = \frac{T}{2}{\text{tr}}\left({{{\mathbf{\Omega }}_{1}}{{\mathbf{\Lambda }}_{1}}{\mathbf{\Lambda }}_{1}^{H}} \right) + \frac{T}{2}{\text{tr}}\left({{{\mathbf{\Omega }}_{2}}{{\mathbf{\Lambda }}_{2}}{\mathbf{\Lambda }}_{2}^{H}} \right){\text{ + }}C, \end{array} $$
(46)

where

$$\begin{array}{*{20}l} {}{{\mathbf{\Omega }}_{1}} = {\left({{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}{{\mathbf{\Sigma }}_{1}}} \right)^{H}}\left({{\alpha_{1}}{{\mathbf{\Lambda }}_{{g_{1}}}} + {\alpha_{2}}{{\mathbf{\Lambda }}_{{g_{2}}}}} \right)\left({{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}{{\mathbf{\Sigma }}_{1}}} \right) - {{\mathbf{I}}_{N}}, \end{array} $$
(47)
$$\begin{array}{*{20}l} {}{{\mathbf{\Omega }}_{2}} = {\left({{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}{{\mathbf{\Sigma }}_{2}}} \right)^{H}}\left({{\alpha_{1}}{{\mathbf{\Lambda }}_{{g_{1}}}} + {\alpha_{2}}{{\mathbf{\Lambda }}_{{g_{2}}}}} \right)\left({{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}{{\mathbf{\Sigma }}_{2}}} \right) - {{\mathbf{I}}_{N}} \end{array} $$
(48)

and

$$\begin{array}{*{20}l} {} C &= \frac{T}{2}{\alpha_{1}}{\text{tr}}\left({{{\mathbf{\Lambda }}_{{g_{1}}}}{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}{{\mathbf{\Lambda }}_{h}}{\mathbf{\Lambda }}_{R}^{H}{\mathbf{\Sigma }}_{g}^{H} + \sigma_{1}^{2}{{\mathbf{I}}_{N}}} \right) \hfill \\ &\quad + \frac{T}{2}{\alpha_{2}}{\text{tr}}\left({{{\mathbf{\Lambda }}_{{g_{2}}}}{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}{{\mathbf{\Lambda }}_{h}}{\mathbf{\Lambda }}_{R}^{H}{\mathbf{\Sigma }}_{g}^{H} + \sigma_{2}^{2}{{\mathbf{I}}_{N}}} \right)\\ &\quad+ E_{1}^{INI} + E_{2}^{INI} \hfill. \end{array} $$
(49)

Define \(\phantom {\dot {i}\!}{\lambda _{i,n}}, {\lambda _{R,n}}, {\lambda _{{h_{i}},n}}, {\lambda _{g,n}}, {\lambda _{{\Psi _{h}},n}}\), and \(\phantom {\dot {i}\!}{\lambda _{{\Psi _{{g_{i}}}},n}}\) as the nth diagonal entries of \(\phantom {\dot {i}\!}{{\mathbf {\Lambda }}_{i}}, {{\mathbf {\Lambda }}_{R}}, {{\mathbf {\Sigma }}_{i}}, {{\mathbf {\Sigma }}_{g}}, {{\mathbf {\Lambda }}_{h}}\), and \(\phantom {\dot {i}\!}{{\mathbf {\Lambda }}_{{g_{i}}}}\), respectively. Following this definitions, the sum retained energy can be transformed in a scalar form as

$$\begin{array}{*{20}l} {}\overset{{\smash{\scriptscriptstyle\frown}}}{E}_{1}^{R} + \overset{{\smash{\scriptscriptstyle\frown}}}{E}_{2}^{R} = \frac{T}{2}\left({\sum\limits_{n = 1}^{N} {{\lambda_{{\Omega_{1}},n}}\lambda_{1,n}^{2}} + \sum\limits_{n = 1}^{N} {{\lambda_{{\Omega_{2}},n}}\lambda_{2,n}^{2}}} \right) + C, \end{array} $$
(50)

where \({\lambda _{{\Omega _{i}},n}}\) denotes the nth diagonal entries of Ωi. With (45) and (50), the original optimization problem (9) with matrix variables is reformulated as the following power allocation problem with scalar variables

$$\begin{array}{*{20}l} &{}\underset{\left({{{\mathbf{\lambda }}_{1}},\;{{\mathbf{\lambda }}_{2}}} \right)}{\max} \;\;\;\frac{T}{2}\left({\sum\limits_{n = 1}^{N} {{\lambda_{{\Omega_{1}},n}}\lambda_{1,n}^{2}} + \sum\limits_{n = 1}^{N} {{\lambda_{{\Omega_{2}},n}}\lambda_{2,n}^{2}}} \right) + C \hfill \end{array} $$
(51a)
$$\begin{array}{*{20}l} & {}s.t.\;\frac{{\sum\limits_{n = 1}^{N} {\lambda_{g,n}^{2}\lambda_{R,n}^{2}\lambda_{{h_{3 - i}},n}^{2}\lambda_{3 - i,n}^{2}} }}{{\sum\limits_{n = 1}^{N} {\lambda_{g,n}^{2}\lambda_{R,n}^{2}{\lambda_{{\Psi_{h}},n}}\sigma_{R}^{2}} {\text{ + }}\sum\limits_{n = 1}^{N} {\frac{{\sigma_{i}^{2}}}{{{\lambda_{{\Psi_{{g_{i}}}},n}}}}} {\text{ + }}\sum\limits_{n = 1}^{N} {\frac{{\sigma_{i,z}^{2}}}{{{\lambda_{{\Psi_{{g_{i}}}},n}}{\rho_{i}}}}} }} \geqslant {\gamma_{i}},i \in \left\{ {1,2} \right\}, \hfill \end{array} $$
(51b)
$$\begin{array}{*{20}l} &{}\;\;\;\;\;\;\sum\limits_{n = 1}^{N} {\lambda_{R,n}^{2}\left({\lambda_{{h_{1}},n}^{2}\lambda_{1,n}^{2}{\text{ + }}\lambda_{{h_{2}},n}^{2}\lambda_{2,n}^{2}{\text{ + }}{\lambda_{{\Psi_{h}},n}}\sigma_{R}^{2}} \right)} \leqslant {P_{R,\;Max}} \hfill, \end{array} $$
(51c)
$$\begin{array}{*{20}l} &{} \;\;\;\;\;\;\sum\limits_{n = 1}^{N} {\lambda_{i,n}^{2}} \leqslant {P_{i,\;Max}},i \in \left\{ {1,2} \right\} \hfill, \end{array} $$
(51d)

where \({{\mathbf {\lambda }}_{i}} = {\left [ {\lambda _{i,1}^{2}, \cdots,\lambda _{i,N}^{2}} \right ]^{T}}\) indicates the allocation power of the subchannels at UE i. Since the linearity of the objective function with the terms of λ1 and λ2, the optimization problem (51) is convex and can be readily solved through the classical optimization scheme, i.e., CVX tool.

5.3 Optimize UAV relay power allocation

Similarly, for given Λ1,Λ2,ρ1, and ρ2, we have the sum retained energy function with UAV relay power allocation matrix ΛR as

$$\begin{array}{*{20}l} \overset{{\smash{\scriptscriptstyle\frown}}}{E}_{1}^{R} + \overset{{\smash{\scriptscriptstyle\frown}}}{E}_{2}^{R} = \frac{T}{2}{\text{tr}}\left({\left({{{\mathbf{\Omega }}_{3}} + {{\mathbf{\Omega }}_{4}}} \right){{\mathbf{\Lambda }}_{R}}{\mathbf{\Lambda }}_{R}^{H}} \right) + \tilde{C}, \end{array} $$
(52)

where

$$\begin{array}{*{20}l} {{\mathbf{\Omega }}_{3}} = {\alpha_{1}}{\mathbf{\Sigma }}_{g}^{H}{{\mathbf{\Lambda }}_{{g_{1}}}}{{\mathbf{\Sigma }}_{g}}\left({{\mathbf{\Pi }} + \sigma_{R}^{2}{{\mathbf{\Lambda }}_{h}}} \right), \end{array} $$
(53)
$$\begin{array}{*{20}l} {{\mathbf{\Omega }}_{4}} = {\alpha_{2}}{\mathbf{\Sigma }}_{g}^{H}{{\mathbf{\Lambda }}_{{g_{2}}}}{{\mathbf{\Sigma }}_{g}}\left({{\mathbf{\Pi }} + \sigma_{R}^{2}{{\mathbf{\Lambda }}_{h}}} \right) \end{array} $$
(54)

and

$$\begin{array}{*{20}l} {}\tilde{C} &= \frac{T}{2}{\text{tr}}\left({{\alpha_{1}}\sigma_{1}^{2}{{\mathbf{I}}_{N}} + {\alpha_{2}}\sigma_{2}^{2}{{\mathbf{I}}_{N}}} \right) + E_{1}^{INI} + E_{2}^{INI}\\ &\quad - \frac{T}{2}{\text{tr}}\left({{{\mathbf{\Lambda }}_{1}}{\mathbf{\Lambda }}_{1}^{H}{\text{ + }}{{\mathbf{\Lambda }}_{2}}{\mathbf{\Lambda }}_{2}^{H}} \right). \end{array} $$
(55)

According to (50), (52) can be transformed in a scalar form as

$$\begin{array}{*{20}l} \overset{{\smash{\scriptscriptstyle\frown}}}{E}_{1}^{R} + \overset{{\smash{\scriptscriptstyle\frown}}}{E}_{2}^{R} = \frac{T}{2}\sum\limits_{n = 1}^{N} {\left({\left({{\lambda_{{\Omega_{3}},n}} + {\lambda_{{\Omega_{4}},n}}} \right)\lambda_{R,n}^{2}} \right)} + \tilde{C}, \end{array} $$
(56)

where \({\lambda _{{\Omega _{3}},n}}\) and \({\lambda _{{\Omega _{4}},n}}\) strand for the nth diagonal entries of Ω3 and Ω4, respectively. Combining (45) with (56), the original optimization problem in (9) can be restated with scalar variables as

$$\begin{array}{*{20}l} &\underset{\left({\;{{\mathbf{\lambda }}_{R}}} \right)} {\max} \;\;\;\frac{T}{2}\sum\limits_{n = 1}^{N} {\left({\left({{\lambda_{{\Omega_{3}},n}} + {\lambda_{{\Omega_{4}},n}}} \right)\lambda_{R,n}^{2}} \right)} + \tilde{C} \hfill \end{array} $$
(57a)
$$\begin{array}{*{20}l} & s.t.\;\sum\limits_{n = 1}^{N} {\left({\left({\lambda_{g,n}^{2}\lambda_{{h_{3 - i}},n}^{2}\lambda_{3 - i,n}^{2} - {\gamma_{i}}\lambda_{g,n}^{2}{\lambda_{{\Psi_{h}},n}}\sigma_{R}^{2}} \right)\lambda_{R,n}^{2}} \right)} \hfill \\ &\; \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\geqslant {\gamma_{i}}\left({\sum\limits_{n = 1}^{N} {\frac{{\sigma_{i}^{2}}}{{{\lambda_{{\Psi_{{g_{i}}}},n}}}} + \sum\limits_{n = 1}^{N} {\frac{{\sigma_{i,z}^{2}}}{{{\lambda_{{\Psi_{{g_{i}}}},n}}{\rho_{i}}}}}} } \right),i \in \left\{ {1,2} \right\}, \hfill \end{array} $$
(57b)
$$\begin{array}{*{20}l} &\;\;\;\;\;\;\sum\limits_{n = 1}^{N} {\left({\left({\lambda_{{h_{1}},n}^{2}\lambda_{1,n}^{2}{\text{ + }}\lambda_{{h_{2}},n}^{2}\lambda_{2,n}^{2}{\text{ + }}{\lambda_{{\Psi_{h}},n}}\sigma_{R}^{2}} \right)\lambda_{R,n}^{2}} \right)} \leqslant {P_{R,\;Max}} \hfill, \end{array} $$
(57c)

where \({{\mathbf {\lambda }}_{R}} = {\left [ {\lambda _{R,1}^{2}, \cdots,\lambda _{R,N}^{2}} \right ]^{T}}\) is the allocation power of the subchannels at UAV relay. Due to the linearity of the objective function with the term of λR, the problem (57) is also convex and can be efficiently solved by CVX tool.

5.4 Optimize PS factors

For given Λ1,Λ2, and ΛR, the objective function of maximizing the sum retained energy with ρ1 and ρ2 can be written as

$$\begin{array}{*{20}l} \overset{{\smash{\scriptscriptstyle\frown}}}{E}_{1}^{R} + \overset{{\smash{\scriptscriptstyle\frown}}}{E}_{2}^{R} = {\hat{\beta_{1}}}\left({1 - {\rho_{1}}} \right){\text{ + }}{\hat{\beta_{2}}}\left({1 - {\rho_{2}}} \right){\text{ + }}{\hat {C}}, \end{array} $$
(58)

where

$$\begin{array}{*{20}l} {}{\hat{\beta_{1}}} = \frac{T}{2}\eta {\text{tr}}\left({{{\mathbf{\Lambda }}_{{g_{1}}}}{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}\left({{\mathbf{\Pi }} + {{\mathbf{\Lambda }}_{h}}\sigma_{R}^{2}} \right){\mathbf{\Lambda }}_{R}^{H}{\mathbf{\Sigma }}_{g}^{H} + \sigma_{1}^{2}{{\mathbf{I}}_{N}}} \right), \end{array} $$
(59)
$$\begin{array}{*{20}l} {}{\hat{\beta_{2}}} = \frac{T}{2}\eta {\text{tr}}\left({{{\mathbf{\Lambda }}_{{g_{2}}}}{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}\left({{\mathbf{\Pi }} + {{\mathbf{\Lambda }}_{h}}\sigma_{R}^{2}} \right){\mathbf{\Lambda }}_{R}^{H}{\mathbf{\Sigma }}_{g}^{H} + \sigma_{2}^{2}{{\mathbf{I}}_{N}}} \right) \end{array} $$
(60)

and

$$\begin{array}{*{20}l} {\hat{C}} = E_{1}^{INI} + E_{2}^{INI} - \frac{T}{2}{\text{tr}}\left({{{\mathbf{\Lambda }}_{1}}{\mathbf{\Lambda }}_{1}^{H}{\text{ + }}{{\mathbf{\Lambda }}_{2}}{\mathbf{\Lambda }}_{2}^{H}} \right). \end{array} $$
(61)

Based on the SNR constraint in (45), we have

$$\begin{array}{*{20}l} {\rho_{i}} \geqslant \frac{{{\gamma_{i}}{\text{tr}}\left({{\mathbf{\Lambda }}_{{g_{i}}}^{- 1}} \right)\sigma_{i,z}^{2}}}{{{{\hat{\mu}}_{i}}}}, \end{array} $$
(62)

where

$$\begin{array}{*{20}l} {} {{\hat{\mu}}_{i}}&{\text{ = tr}}\left({{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}{{\mathbf{\Sigma }}_{3 - i}}{{\mathbf{\Lambda }}_{3 - i}}{\mathbf{\Lambda }}_{3 - i}^{H}{\mathbf{\Sigma }}_{3 - i}^{H}{\mathbf{\Lambda }}_{R}^{H}{\mathbf{\Sigma }}_{g}^{H}} \right) \hfill \\ & \quad- {\gamma_{i}}\left({{\text{tr}}\left({{\mathbf{\Lambda }}_{R}^{H}{\mathbf{\Sigma }}_{g}^{H}{{\mathbf{\Sigma }}_{g}}{{\mathbf{\Lambda }}_{R}}\;{{\mathbf{\Lambda }}_{h}}} \right)\sigma_{R}^{2}{\text{ + tr}}\left({{\mathbf{\Lambda }}_{{g_{i}}}^{- 1}} \right)\sigma_{i}^{2}} \right) \hfill. \end{array} $$
(63)

Combining (58) and (62), the original optimization problem can be reformulated as

$$\begin{array}{*{20}l} &\underset{\left({{\rho_{1}},\;{\rho_{2}}} \right)}{\max} \;\;\;{{\hat{\beta}}_{1}}\left({1 - {\rho_{1}}} \right){\text{ + }}{{\hat{\beta}}_{2}}\left({1 - {\rho_{2}}} \right){\text{ + }}{\hat{C}} \hfill \end{array} $$
(64a)
$$\begin{array}{*{20}l} &s.t.\;\;{\rho_{i}} \geqslant \frac{{{\gamma_{i}}{\text{tr}}\left({{\mathbf{\Lambda }}_{{g_{i}}}^{- 1}} \right)\sigma_{i,z}^{2}}}{{{{\hat{\mu}}_{i}}}},i \in \left\{ {1,2} \right\}, \hfill \end{array} $$
(64b)
$$\begin{array}{*{20}l} &\;\;\;\;\;\;\;0 < {\rho_{i}} < 1\;,i \in \left\{ {1,2} \right\}. \hfill \end{array} $$
(64c)

From the relationship of variates in (64), the optimal closed-form expression of PS factors can be expressed as

$$\begin{array}{*{20}l} {\rho_{i,\;opt}} = \frac{{{\gamma_{i}}{\text{tr}}\left({{\mathbf{\Lambda }}_{{g_{i}}}^{- 1}} \right)\sigma_{i,z}^{2}}}{{{{\hat{\mu}}_{i}}}},i \in \left\{ {1,2} \right\}. \end{array} $$
(65)

In summary, the proposed low-complexity AO scheme based on GSVD for maximizing the sum remained energy can be summarized as Algorithm 2.

Complexity Analysis: Note that the problem (45) can be effectively solved by three convex problems with scalar variables (51), (57), and (64). Based on the computational cost analysis in [18], the complexity of the low-complexity beamforming design for maximizing the sum retained energy is about \(\mathcal {O}\left ({{M_{Iter}} \cdot {N^{3.5}}} \right)\). It is obvious that Algorithm 2 has lower computational complexity than Algorithm 1. Therefore, Algorithm 2 is more suitable for large-scale antennas system.

6 Results and discussion

In this section, some simulation results are provided to evaluate the effectiveness of the proposed schemes for UAV-enabled two-way relay system with SWIPT. As illustrated in Fig. 2, the UAV relay hovers at a constant altitude that is set as dv=5 m and the horizontal distance between UE 1 and UE 2 is fixed as dh=30 m. Moreover, we denote the horizontal distances of UE 1-to-relay and relay-to- UE 2 as ξdh and (1−ξ)dh, respectively. ξ∈(0, 1) is the location ratio of UAV relay. It is intuitive that the real distances of UE 1-to-relay and relay-to- UE 2 can be indicated as \({L_{1}} = \sqrt {d_{v}^{2} + {{\left ({\xi {d_{h}}} \right)}^{2}}}\) and \({L_{2}} = \sqrt {d_{v}^{2} + {{\left ({\left ({1 - \xi } \right){d_{h}}} \right)}^{2}}} \), respectively. Following [21, 30], we model the MIMO channel matrixes as \({{\mathbf {H}}_{i}} = {\mathbf {\breve{H} }}_{i}\sqrt {L_{i}^{- \varepsilon }} \) and \({{\mathbf {G}}_{i}} = {\mathbf {\breve{G}}}_{i}\sqrt {L_{i}^{- \varepsilon }} \), where \(L_{i}^{- \varepsilon }\) strands for the large-scale fading, H ̌i and G ̌i indicate the small-scale Rician fading (Rician factor K=5). Similar to [17], the path loss exponent of the high altitude transmission is set as ε=2. For all simulations, we assume that the noise variances at all transceivers are same as \(\sigma _{R}^{2}{\text { = }}\sigma _{i,z}^{2}{\text { = }}\sigma _{i}^{2}{\text { = }} - 60\) dBm. We also specify η=0.8,ξ=0.5,N=2,T=1 second, ρi=0.5,γi=34 dB, \(E_{i}^{INI} = 50\) mJ, Pi, Max=20 dBm and PR, Max=20 dBm unless otherwise explained.

Fig. 2
figure 2

The location of UAV-enable two-way relay system

Figure 3 reveals the convergence performance of the proposed AO scheme and low-complexity GSVD scheme. It is found that the curve of the proposed AO scheme is similar to that of low-complexity GSVD scheme. They take no more than 5 iterations to converge. Compared with the AO scheme, the low-complexity GSVD scheme transforms the original optimization problem with matrix variables into the power allocation optimization problem with scalar variables. However, two schemes have a similar iteration process. Interestingly, the proposed AO scheme has a better performance of the sum retained energy than that of low-complexity GSVD scheme.

Fig. 3
figure 3

Sum retained energy versus the number of iterations

Figure 4 depicts the various performance results of different schemes versus the maximum transmission power of UAV relay. It is intuitive from Fig. 4a that both the proposed AO scheme and low-complexity GSVD scheme significantly outperform the average power assignment (APA) scheme in [31]. This is because compared to APA scheme, the proposed AO scheme and low-complexity GSVD scheme can take full advantage of the robust beamforming matrixes and the effective PS factors to match the MIMO channels. Furthermore, the performance gaps between two proposed schemes and APA scheme grow as the maximum transmission power of UAV relay increases. Besides, it is found that the gap between the proposed AO scheme and low-complexity GSVD scheme is quite small due to the similar iterative process.

Fig. 4
figure 4

Various performance results of different schemes versus the maximum transmission power of UAV relay (dBm) PR,Max. a Sum remained energy versus PR,Max; b Sum cost energy versus PR,Max; c Sum harvested energy versus PR,Max

In Fig. 4b, we check the influence of the maximum transmission power of UAV relay on the sum cost energy of three schemes. It is seen that the APA scheme needs to cost more energy for the same SNR threshold than the proposed AO scheme and low-complexity GSVD scheme at low PR,Max. Moreover, as expected, the cost energy of APA scheme decreases as PR,Max increases. This is because for the fixed SNR requirement and given PS factors, APA scheme requires less transmission power of UE to guarantee the SNR threshold as PR,Max grows.

In Fig. 4c, we investigate the impact of PR,Max on the sum harvested energy of three different schemes. It is observed that the proposed AO scheme and low-complexity GSVD scheme harvest more energy than APA scheme at high PR,Max. This happens because compared to APA scheme, two proposed schemes transfer more energy in MAC phase which results in the higher received power at UEs in BC phase. Besides, two proposed schemes have more efficient beamforming designs and more suitable PS factors to satisfy the SNR thresholds. Specifically, the performance gap of harvested energy between the proposed AO scheme and low-complexity GSVD scheme allocation factor algorithm is small due to the similar iteration procedure.

Figure 5 illustrates the influence of the location ratio of UAV relay ξ on the sum remained energy of three schemes. Obviously, the proposed AO scheme and low-complexity GSVD scheme have better performance than APA scheme. In particular, the better remained energy performance of two proposed schemes are obtained when the UAV relay is closer to UEs. One reason is that the closer UAV relay is located to UEs, the more EH efficiency UEs have. In addition, when UAV relay is closer to UE, the MIMO channel between UAV relay and UE has better channel quality. This is consistent with the simulation result described in [18].

Fig. 5
figure 5

Sum remained energy versus the location ratio of UAV relay ξ,PR,Max=25 dBm

7 Conclusion

In this paper, the joint transceiver beamforming and PS factor optimization of UAV-enabled two-way relaying system have been investigated. SWIPT technique was considered at UEs to prolong the battery life. We proposed the AO scheme and the low-complexity GSVD scheme to maximize the sum remained energy of two UEs. It was found that both the proposed AO scheme and low-complexity GSVD scheme significantly outperform APA scheme with fixed PS factors. In particular, the low-complexity GSVD scheme provides an outstanding balance between complexity and performance. More importantly, it is verified that the low-complexity GSVD scheme performs close to the proposed AO scheme. Besides, we conclude that the location ratio of UAV relay has a huge influence on the sum remained energy.