%% BioMed_Central_Tex_Template_v1.06
%%                                      %
%  bmc_article.tex            ver: 1.06 %
%                                       %

%%IMPORTANT: do not delete the first line of this template
%%It must be present to enable the BMC Submission system to
%%recognise this template!!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                     %%
%%  LaTeX template for BioMed Central  %%
%%     journal article submissions     %%
%%                                     %%
%%          <8 June 2012>              %%
%%                                     %%
%%                                     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                 %%
%% For instructions on how to fill out this Tex template           %%
%% document please refer to Readme.html and the instructions for   %%
%% authors page on the biomed central website                      %%
%% http://www.biomedcentral.com/info/authors/                      %%
%%                                                                 %%
%% Please do not use \input{...} to include other tex files.       %%
%% Submit your LaTeX manuscript as one .tex document.              %%
%%                                                                 %%
%% All additional figures and files should be attached             %%
%% separately and not embedded in the \TeX\ document itself.       %%
%%                                                                 %%
%% BioMed Central currently use the MikTex distribution of         %%
%% TeX for Windows) of TeX and LaTeX.  This is available from      %%
%% http://www.miktex.org                                           %%
%%                                                                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% additional documentclass options:
%  [doublespacing]
%  [linenumbers]   - put the line numbers on margins

%%% loading packages, author definitions

%\documentclass[twocolumn]{bmcart}% uncomment this for twocolumn layout and comment line below
\documentclass{bmcart}

%%% Load packages
%\usepackage{amsthm,amsmath}
%\RequirePackage{natbib}
%\RequirePackage[authoryear]{natbib}% uncomment this for author-year bibliography
%\RequirePackage{hyperref}
\usepackage[utf8]{inputenc} %unicode support
\usepackage{amsmath,amssymb,amsfonts}
\usepackage{algorithmic}
\usepackage{graphicx}
\usepackage{textcomp}
\usepackage{color}
\usepackage{subcaption}
\usepackage{algorithm}
\usepackage{algorithmic}

\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{bm}
\usepackage{color}
\usepackage{stfloats}
\usepackage{multicol}
\usepackage{CJK}
\usepackage{epstopdf}

\usepackage{algorithm}
\usepackage{algorithmic}

%\renewcommand{\algorithmicrequire}{\textbf{Input:}} % Use Input in the format of Algorithm
%\renewcommand{\algorithmicensure}{\textbf{Output:}} % Use Output in the format of Algorithm
%\usepackage[applemac]{inputenc} %applemac support if unicode package fails
%\usepackage[latin1]{inputenc} %UNIX support if unicode package fails


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                             %%
%%  If you wish to display your graphics for   %%
%%  your own use using includegraphic or       %%
%%  includegraphics, then comment out the      %%
%%  following two lines of code.               %%
%%  NB: These line *must* be included when     %%
%%  submitting to BMC.                         %%
%%  All figure files must be submitted as      %%
%%  separate graphics through the BMC          %%
%%  submission process, not included in the    %%
%%  submitted article.                         %%
%%                                             %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\def\includegraphic{}
\def\includegraphics{}



%%% Put your definitions there:
\startlocaldefs
\endlocaldefs


%%% Begin ...
\begin{document}

%%% Start of article front matter
\begin{frontmatter}

\begin{fmbox}
\dochead{Research}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the title of your article here     %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{Design and Optimization for UAV-Enabled Two-Way Relaying System with SWIPT}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors here                   %%
%%                                          %%
%% Specify information, if available,       %%
%% in the form:                             %%
%%   <key>={<id1>,<id2>}                    %%
%%   <key>=                                 %%
%% Comment or delete the keys which are     %%
%% not used. Repeat \author command as much %%
%% as required.                             %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\author[
   addressref={aff1},
   email={wjl200609@gmail.com}
]{\inits{J}\fnm{Jinlong} \snm{Wang}}
\author[
   addressref={aff1},                   % id's of addresses, e.g. {aff1,aff2}
   corref={aff1},                       % id of corresponding address, if any
%   noteref={n1},                        % id's of article notes, if any
   email={gwang51@hit.edu.cn}   % email address
]{\inits{G}\fnm{Gang} \snm{Wang}}
\author[
   addressref={aff1},
   email={chenguanyi@gmail.com}
]{\inits{G}\fnm{Guanyi} \snm{Chen}}
\author[
   addressref={aff1},
   email={libo1983@hit.edu.cn}
]{\inits{B}\fnm{Bo} \snm{Li}}
\author[
   addressref={aff1},
   email={liuchungang@hebtu.edu.cn}
]{\inits{C}\fnm{Ruofei} \snm{Zhou}}
\author[
   addressref={aff2},
   email={ryzhang19@njust.edu.cn}
]{\inits{R}\fnm{Ruoyu} \snm{Zhang}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors' addresses here        %%
%%                                          %%
%% Repeat \address commands as much as      %%
%% required.                                %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\address[id=aff1]{%                           % unique id
  \orgname{Communication Research Center, Harbin Institute of Technology}, % university, etc
  %\street{Waterloo Road},                     %
  \postcode{150001}                                % post or zip code
  \city{Harbin},                              % city
  \cny{China}                                    % country
}

\address[id=aff2]{%                           % unique id
  \orgname{Ministerial Key Laboratory of JGMT, Nanjing University of Science and Technology}, % university, etc
  %\street{Waterloo Road},                     %
  \postcode{210094}                                % post or zip code
  \city{Nanjing},                              % city
  \cny{China}                                    % country
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter short notes here                   %%
%%                                          %%
%% Short notes will be after addresses      %%
%% on first page.                           %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\begin{artnotes}
%%\note{Sample of title note}     % note to the article
%\note[id=n1]{Equal contributor} % note, connected to author
%\end{artnotes}

\end{fmbox}% comment this for two column layout

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Abstract begins here                 %%
%%                                          %%
%% Please refer to the Instructions for     %%
%% authors on http://www.biomedcentral.com  %%
%% and include the section headings         %%
%% accordingly for your article type.       %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{abstractbox}

\begin{abstract} % abstract
In this paper, we investigate the resource allocation scheme for an unmanned-aerial-vehicle-enable (UAV-enabled) two-way relaying system with simultaneous wireless information and power transfer (SWIPT), where two user equipments (UEs) exchange information with the help of UAV relay and harvest energy through power splitting (PS) scheme. Under the transmission power constraints at UEs and UAV relay, a non-convex intractable optimization problem is formulated which maximizes the sum retained energy of two UEs while satisfying the minimum signal-to-noise ratio requirement. We decouple the complicated beamforming and PS factors optimization problem into three solvable subproblems and propose an efficient alternating optimization scheme. Subsequently, in order to reduce the complexity, a robust scheme based on generalized singular value decomposition (GSVD) is designed. Finally, numerical results verify the robustness and effectiveness of two proposed schemes.
\end{abstract}
%\parttitle{First part title} %if any
%Text for this section.
%
%\parttitle{Second part title} %if any
%Text for this section.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The keywords begin here                  %%
%%                                          %%
%% Put each keyword in separate \kwd{}.     %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{keyword}
\kwd{UAV-enabled}
\kwd{resource allocation}
\kwd{SWIPT}
\kwd{MIMO}
\kwd{GSVD}
\end{keyword}

% MSC classifications codes, if any
%\begin{keyword}[class=AMS]
%\kwd[Primary ]{}
%\kwd{}
%\kwd[; secondary ]{}
%\end{keyword}

\end{abstractbox}
%
%\end{fmbox}% uncomment this for twcolumn layout

\end{frontmatter}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Main Body begins here                %%
%%                                          %%
%% Please refer to the instructions for     %%
%% authors on:                              %%
%% http://www.biomedcentral.com/info/authors%%
%% and include the section headings         %%
%% accordingly for your article type.       %%
%%                                          %%
%% See the Results and Discussion section   %%
%% for details on how to create sub-sections%%
%%                                          %%
%% use \cite{...} to cite references        %%
%%  \cite{koon} and                         %%
%%  \cite{oreg,khar,zvai,xjon,schn,pond}    %%
%%  \nocite{smith,marg,hunn,advi,koha,mouse}%%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%
%% Background %%
%%
%\section*{Content}
%Text and results for this section, as per the individual journal's instructions for authors. %\cite{koon,oreg,khar,zvai,xjon,schn,pond,smith,marg,hunn,advi,koha,mouse}


%\nocite{oreg,schn,pond,smith,marg,hunn,advi,koha,mouse}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Backmatter begins here                   %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
Recently, unmanned aerial vehicles (UAVs) have received dramatical attentions in wireless communication networks since they have the capability of flexible mobility [1-4]. Compared to 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 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. [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 equipments (UEs) and forwards the result to the access point. 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 [11-14]. 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 [15-19]. 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 performance [28-30]. In [28], 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. [29] considered the throughput maximization problem in a UAV relay system with time-sharing mechanism, where the UAV relay is capable of SWIPT. In [30], 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.

\section{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.
\begin{itemize}
\item We propose a UAV-enabled two-way relaying system with SWIPT. The system model distinguishes from the existing relaying system, since the UEs are capable of PS protocol and multi-antennas are equipped at all terminals.
\item 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.
\item 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.
\item 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.
\end{itemize}

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

\section{System model and problem formulation}
\subsection{System Model}
As shown in Figure 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., ${\text{UE}}\;1$ and  ${\text{UE}}\;{\text{2}}$, are installed with the same antennas as UAV relay. Due to the barriers or buildings, there is no direct link between ${\text{UE}}\;1$ and ${\text{UE}}\;{\text{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 (AF) way is applied at UAV relay.

With the half-duplex relay, the communication period block $T$ is equally partitioned into two separate phase, 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{align}
{{\mathbf{y}}_{\text{R}}} = {{\mathbf{H}}_{\text{1}}}{{\mathbf{W}}_1}{{\mathbf{x}}_1} + {{\mathbf{H}}_2}{{\mathbf{W}}_2}{{\mathbf{x}}_2} + {{\mathbf{n}}_R},
\label{eq:1}
%\tag{12}
\end{align}
where ${{\mathbf{H}}_i} \in {\mathbb{C}^{N \times N}}$ is the MIMO channel matrix from ${\text{UE}}\;i$ to UAV relay node, ${{\mathbf{W}}_i} \in {\mathbb{C}^{N \times N}}$ represents the beamforming matrix of ${\text{UE}}\;i$, which should be under the transmission power constraints ${\text{tr}}\left( {{{\mathbf{W}}_i}{{\mathbf{W}}_i}^H} \right) = {P_i}$ and ${P_i} \leqslant {P_{i,\;Max}}$, where ${P_i}$ and ${P_{i,\;Max}}$ denote the realistic transmission power and the limitation of maximum transmission power at ${\text{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 ${\text{tr}}\left( {{{\mathbf{W}}_R}{\mathbf{R}}{{\mathbf{W}}_R}^H} \right) = {P_R}$ and ${P_R} \leqslant {P_{R,\;Max}}$, where ${\mathbf{R}} = {{\mathbf{H}}_{\text{1}}}{{\mathbf{W}}_{\text{1}}}{\mathbf{W}}_1^H{\mathbf{H}}_1^H + {{\mathbf{H}}_{\text{2}}}{{\mathbf{W}}_{\text{2}}}{\mathbf{W}}_2^H{\mathbf{H}}_2^H + \sigma _R^2{{\mathbf{I}}_N}$, ${P_R}$ and ${P_{R,\;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{align}
{{\mathbf{x}}_R} = {{\mathbf{W}}_R}{{\mathbf{y}}_R} = {{\mathbf{W}}_R}\left( {{{\mathbf{H}}_{\text{1}}}{{\mathbf{W}}_{\text{1}}}{{\mathbf{x}}_{\text{1}}} + {{\mathbf{H}}_{\text{2}}}{{\mathbf{W}}_{\text{2}}}{{\mathbf{x}}_{\text{2}}} + {{\mathbf{n}}_R}} \right).
\label{eq:2}
%\tag{12}
\end{align}

The received RF signal at ${\text{UE}}\;i$ is split into two portions with a PS factor, i.e., ${\rho _i} \in \left( {0,1} \right)$. ${\rho _i}$ portion is used for ID, and $1 - \rho _i $ portion is for EH. In order to reduce the complexity of MIMO relay system, we assume that PS factors of all antennas at ${\text{UE}}\;i$ are same. The ID signal at ${\text{UE}}\;i$ can be written as
\begin{align}
{{\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}} + {{\mathbf{G}}_i}{{\mathbf{W}}_R}{{\mathbf{n}}_R} + {{\mathbf{n}}_i}} \right) + {{\mathbf{n}}_{i,z}},
\label{eq:3}
\end{align}
where ${{\mathbf{n}}_i}$ indicates the zero-mean AWGN vector at ${\text{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 ${{\mathbf{G}}_i}{{\mathbf{W}}_R}{{\mathbf{H}}_i}{{\mathbf{W}}_i}{{\mathbf{x}}_i}$ can be integrally canceled. Therefore, the received signal at ${\text{UE}}\;i$ can be further obtained as
\begin{align}
{{\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}}.
\label{eq:4}
\end{align}

According to (4), the received SNR at ${\text{UE}}\;i$ can be represented by
\begin{align}
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}}}.
\label{eq:5}
\end{align}

The harvested energy at ${\text{UE}}\;i$ can be given by
\begin{align}
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),
\label{eq:6}
\end{align}
where $\eta  \in \left( {0,1} \right)$ is  the energy conversion efficiency of ${\text{UE}}\;i$. It is worth pointing out that the self-interference term ${{\mathbf{G}}_i}{{\mathbf{W}}_R}{{\mathbf{H}}_i}{{\mathbf{W}}_i}{{\mathbf{x}}_i}$  can also be utilized for EH.

\subsection{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 ${\text{UE}}\;i$ can be expressed as
\begin{align}
E_i^R = E_i^{INI}{\text{ + }}E_{i,e}^H - \frac{T}{2}{P_i},
\label{eq:7}
\end{align}
where $E_i^{INI}$ denotes the initial energy at ${\text{UE}}\;i$. The sum retained energy of two UEs can be expressed as
\begin{align}
E_{sum}^R = \sum\limits_{i = 1}^2 {E_i^R}.
\label{eq:8}
\end{align}

Based on (5) and (8), the optimization problem of sum retained energy is formulated as
\begin{align}
 &\mathop {\max }\limits_{\left( {{{\mathbf{F}}_1},\;{{\mathbf{F}}_2},\;{{\mathbf{F}}_R},\;{\rho _1},\;{\rho _2}} \right)} \;\;\;E_{sum}^R \hfill
  \tag{9a}
  \\
  &s.t.\;\;SN{R_i} \geqslant {\gamma _i},i \in \left\{ {1,2} \right\} \hfill,
 \tag{9b}
  \\
  &\;\;\;\;\;\;\operatorname{Tr} \left( {{{\mathbf{W}}_R}{\mathbf{RW}}_R^H} \right) \leqslant {P_{R,\;Max}} \hfill,
   \tag{9c}
  \\
  &\;\;\;\;\;\;\operatorname{Tr} \left( {{{\mathbf{W}}_i}{{\mathbf{W}}_i}^H} \right) \leqslant {P_{i,\;Max}},i \in \left\{ {1,2} \right\} \hfill,
   \tag{9d}
  \\
  &\;\;\;\;\;\;\;0 < {\rho _i} < 1\;,i \in \left\{ {1,2} \right\} \hfill ,
  \tag{9e}
\end{align}
where ${\gamma _i}$ indicates SNR threshold of ${\text{UE}}\;i$, constraint (9b) represents the minimum SNR requirement of ${\text{UE}}\;i$,  (9c) and (9d) denote the maximum transmission power of UAV relay and ${\text{UE}}\;i$, respectively, and (9e) is the PS factor constraint of ${\text{UE}}\;i$. Since the optimization problem is non-convex, it is intractable to obtain the optimal solution of maximizing the sum retained energy.

\section{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 UEs beamforming matrixes ${{\mathbf{W}}_1}$ and ${{\mathbf{W}}_2}$ with fixed ${{\mathbf{W}}_R}$, ${\rho _1}$ and ${\rho _2}$. Secondly, we try to obtain the optimal relay beamforming matrix ${{\mathbf{W}}_R}$ by assuming the other parameters are given. Finally, the optimal PS factors, i.e., ${\rho _1}$ and ${\rho _2}$, are investigated.

\subsection{Optimize two UEs beamforming matrixes (${{\mathbf{W}}_1}$ and ${{\mathbf{W}}_2}$)}
With the above analysis, we first design the optimal solution of ${{\mathbf{W}}_1}$ and ${{\mathbf{W}}_2}$ for maximizing the sum retained energy. Substituting ${\mathbf{R}} = {{\mathbf{H}}_{\text{1}}}{{\mathbf{W}}_{\text{1}}}{\mathbf{W}}_1^H{{\mathbf{H}}_{\text{1}}}^H + {{\mathbf{H}}_{\text{2}}}{{\mathbf{W}}_{\text{2}}}{\mathbf{W}}_2^H{{\mathbf{H}}_{\text{2}}}^H + \sigma _R^2{{\mathbf{I}}_N}$ into (9), we have the objective function as
\begin{align}
  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 \nonumber\\
  & + \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 \hfill ,
  \tag{10}
\end{align}
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}$. Base on the trace property ${\text{tr}}\left( {{\mathbf{AB}}} \right) = {\text{tr}}\left( {{\mathbf{BA}}} \right)$, the total retained energy can be rewritten as
\begin{align}
  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\nonumber \\
   &+ \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,
  \tag{11}
\end{align}
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{align}
 {\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),
  \tag{12}
\end{align}
where $vec\left(  \cdot  \right)$ is the matrix vectorization operator, $\otimes $ 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{align}
  {\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 \nonumber\\
   &= {\mathbf{w}}_1^H\left( {{{\mathbf{I}}_N} \otimes {{\mathbf{\Psi }}_1}} \right){{\mathbf{w}}_1} \hfill \nonumber\\
  &{\text{ = tr}}\left( {\left( {{{\mathbf{I}}_N} \otimes {{\mathbf{\Psi }}_1}} \right){{\mathbf{w}}_1}{\mathbf{w}}_1^H} \right) \hfill ,
  \tag{13}
\end{align}
and
\begin{align}
{\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),
  \tag{14}
\end{align}
where ${{\mathbf{w}}_1} = vec\left( {{{\mathbf{W}}_1}} \right)$ . 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{align}
  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 \nonumber\\
  & + \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,
  \tag{15}
\end{align}
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{align}
{\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},
  \tag{16}
\end{align}
\begin{align}
{\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},
  \tag{17}
\end{align}
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{align}
{\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}}.
  \tag{18}
\end{align}

Plugging (16), (17) and (18) back into (9), the original optimization problem can be reformulated as
\begin{align}
  &\mathop {\max }\limits_{\left( {{{\mathbf{\Phi }}_1},\;{{\mathbf{\Phi }}_2} \succcurlyeq {\mathbf{0}}} \right)} \;\;\;\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 \hfill
 \tag{19a}
  \\
  &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} \hfill,
  \tag{19b}
  \\
  &\;\;\;\;\;\;\;{\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,
   \tag{19c}
   \\
 & \;\;\;\;\;\;\;{\text{tr}}\left( {{{\mathbf{\Phi }}_i}} \right) \leqslant {P_{i,\;Max}},\;i \in \left\{ {1,2} \right\} \hfill ,\tag{19d}
  \\
 & \;\;\;\;\;\;\;{\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}} \hfill ,
  \tag{19e}
\end{align}
where ${{\mathbf{Q}}_1}{\text{ = }}\left( {{{\mathbf{I}}_N} \otimes {{\mathbf{\Psi }}_1}} \right) + \left( {{{\mathbf{I}}_N} \otimes {{\mathbf{\Psi }}_2}} \right) - {{\mathbf{I}}_E}$ and ${{\mathbf{Q}}_2}{\text{ = }}\left( {{{\mathbf{I}}_N} \otimes {{\mathbf{\Psi }}_3}} \right) + \left( {{{\mathbf{I}}_N} \otimes {{\mathbf{\Psi }}_4}} \right) - {{\mathbf{I}}_E}$. According to [24], the optimization problem (19) can be converted into a semidefinite 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{*}}$.

\subsection{Optimize UAV relay beamforming (${{\mathbf{W}}_R}$)}
Similarly, with given ${{\mathbf{W}}_1}$, ${{\mathbf{W}}_2}$ and two UEs PS factors, i.e., ${\rho _1}$, ${\rho _2}$, the total retained energy of two UEs is equivalent to
\begin{align}
E_{sum}^R = \frac{T}{2}{\text{tr}}\left( {{{{\mathbf{\tilde \Psi }}}_1}{{\mathbf{W}}_R}{\mathbf{RW}}_R^H + {{{\mathbf{\tilde \Psi }}}_2}{{\mathbf{W}}_R}{\mathbf{RW}}_R^H} \right) + \tilde J,
  \tag{20}
\end{align}
where ${{\mathbf{\tilde \Psi }}_1} = {\alpha _1}{\mathbf{G}}_1^H{{\mathbf{G}}_1}$, ${{\mathbf{\tilde \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{align}
E_{sum}^R = \frac{T}{2}{\text{tr}}\left( {\left( {\left( {{{\mathbf{R}}^T} \otimes {{{\mathbf{\tilde \Psi }}}_1}} \right) + \left( {{{\mathbf{R}}^T} \otimes {{{\mathbf{\tilde \Psi }}}_2}} \right)} \right){{\mathbf{w}}_R}{\mathbf{w}}_R^H} \right) + \tilde J,
  \tag{21}
\end{align}
where ${{\mathbf{w}}_R} = vec\left( {{{\mathbf{W}}_R}} \right)$. Further, the constraints in (9) can be rewritten as
\begin{align}
\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),
  \tag{22}
\end{align}
\begin{align}
\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),
  \tag{23}
\end{align}
and
\begin{align}
\text{tr}\left( {\left( {{{\mathbf{R}}^T} \otimes {{\mathbf{I}}_N}} \right){{\mathbf{w}}_R}{\mathbf{w}}_R^H} \right) \leqslant {P_{R,\;Max}},
  \tag{24}
\end{align}
where
\begin{align}
{{\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),
  \tag{25}
\end{align}

\begin{align}
{{\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).
  \tag{26}
\end{align}

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{align}
  &\mathop {\max }\limits_{\left( {{{\mathbf{\Phi }}_R} \succcurlyeq {\mathbf{0}}} \right)} \;\frac{T}{2}{\text{tr}}\left( {\left( {\left( {{{\mathbf{R}}^T} \otimes {{{\mathbf{\tilde \Psi }}}_1}} \right) + \left( {{{\mathbf{R}}^T} \otimes {{{\mathbf{\tilde \Psi }}}_2}} \right)} \right){{\mathbf{\Phi }}_R}} \right) + \tilde J \hfill
  \tag{27a}
  \\
 & 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,
  \tag{27b}
  \\
& \;\;\;\;\;\;\; \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,
    \tag{27c}
  \\
 &\;\;\;\;\;\;\; \text{tr}\left( {\left( {{{\mathbf{R}}^T} \otimes {{\mathbf{I}}_N}} \right){{\mathbf{\Phi }}_R}} \right) \leqslant {P_{R,\;Max}} \hfill.
  \tag{27d}
\end{align}

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{*}}$.

\subsection{Optimize two PS factors (${\rho _1}$, ${\rho _1}$)}
For given ${{\mathbf{W}}_1}$, ${{\mathbf{W}}_2}$ and ${{\mathbf{W}}_R}$, the objective function of maximizing the sum retained energy with ${\rho _1}$ and ${\rho _2}$ can be written as
\begin{align}
E_{sum}^R = {\beta _1}\left( {1 - {\rho _1}} \right){\text{ + }}{\beta _2}\left( {1 - {\rho _2}} \right){\text{ + }}\hat J,
  \tag{28}
\end{align}
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{align}
{\rho _1} \geqslant \frac{{\sigma _{1,z}^2{\gamma _1}}}{{{\mu _1}}},
  \tag{29}
\end{align}
\begin{align}
{\rho _2} \geqslant \frac{{\sigma _{2,z}^2{\gamma _2}}}{{{\mu _2}}},
  \tag{30}
\end{align}
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{align}
  &\mathop {\max }\limits_{\left( {{\rho _1},{\rho _2}} \right)} \;{\beta _1}\left( {1 - {\rho _1}} \right){\text{ + }}{\beta _2}\left( {1 - {\rho _2}} \right){\text{ + }}\hat J \hfill
   \tag{31a}
  \\
 & s.t.\;\;\;{\rho _1} \geqslant \frac{{\sigma _{1,z}^2{\gamma _1}}}{{{\mu _1}}} \hfill,
  \tag{31b}
  \\
  &\;\;\;\;\;\;\;\;{\rho _2} \geqslant \frac{{\sigma _{2,z}^2{\gamma _2}}}{{{\mu _2}}} \hfill,
  \tag{31c}
  \\
 & \;\;\;\;\;\;\;\;0 < {\rho _1} < 1 \hfill,
  \tag{31d}
   \\
 & \;\;\;\;\;\;\;\;0 < {\rho _2} < 1 \hfill.
  \tag{31e}
\end{align}

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{align}
{\rho _{1,\;opt}} = \frac{{\sigma _{1,z}^2{\gamma _1}}}{{{\mu _1}}},
  \tag{32}
\end{align}
\begin{align}
{\rho _{2,\;opt}} = \frac{{\sigma _{2,z}^2{\gamma _2}}}{{{\mu _2}}}.
  \tag{33}
\end{align}

In summary, the proposed AO scheme for maximizing the sum remained energy can be summarized as Algorithm 1.
\begin{algorithm}[htb]
\caption{The proposed AO scheme for maximizing the sum retained energy}
\label{alg:example}
\begin{algorithmic}[1] %这个1 表示每一行都显示数字
%\REQUIRE The transmit power vector ${\bf{p}}$, the step-size $r$, constant\\
%\quad \; $\delta $ and $\varepsilon $. \\ %算法的输入参数：Input
%\ENSURE The optimal transmit power vector ${{\bf{p}}^*}$.  \\ %算法的输出：Output
\STATE {\bf{Initialization:}}  Set ${\rho _1} = {\rho _2} = 0.5$, ${{\mathbf{W}}_R} = \sqrt {{{{P_{R,\;Max}}} \mathord{\left/
 {\vphantom {{{P_{R,\;Max}}} N}} \right.
 \kern-\nulldelimiterspace} N}} {{\mathbf{I}}_N}$, ${{\mathbf{W}}_1} = \sqrt {{{{P_{1,\;Max}}} \mathord{\left/
 {\vphantom {{{P_{1,\;Max}}} N}} \right.
 \kern-\nulldelimiterspace} N}} {{\mathbf{I}}_N}$ and ${{\mathbf{W}}_2} = \sqrt {{{{P_{2,\;Max}}} \mathord{\left/
 {\vphantom {{{P_{2,\;Max}}} N}} \right.
 \kern-\nulldelimiterspace} N}} {{\mathbf{I}}_N}$.  % while循环，需要以\ENDWHILE作为结尾

\STATE {\bf{Repeat}\\
\texttt{a) Update ${{\mathbf{\Phi }}_1}$ and ${{\mathbf{\Phi }}_2}$  by solving the convex problem (19) and to obtain ${{\mathbf{W}}_1}$ and ${{\mathbf{W}}_2}$;\\
b) Update ${{\mathbf{\Phi }}_R}$ by solving the convex problem (27) and to obtain ${{\mathbf{W}}_R}$;\\
c) Update the optimal PS factors solution ${\rho _{1,\;opt}}$ and ${\rho _{2,\;opt}}$.}}
\STATE {\bf{Until convergence}}.
    %return返回值
\end{algorithmic}
\end{algorithm}

\emph{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 ${M_{Iter}}$ indicates the number of iterations.

\section{Low-complexity beamforming design based on GSVD technique}
The proposed AO scheme of UAV-enable two-way relay system in Section 3 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.
\subsection{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 ${{\mathbf{H}}_1}$ and ${{\mathbf{H}}_2}$ as
\begin{align}
{{\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,
  \tag{34}
\end{align}
where ${{\mathbf{B}}_h}$ 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{align}
{{\mathbf{G}}_{BC}} = {{\mathbf{B}}_g}{{\mathbf{\Sigma }}_g}{\mathbf{A}}_g^H,
  \tag{35}
\end{align}
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), ${{\mathbf{G}}_1}$ and ${{\mathbf{G}}_2}$ are decomposed as [17, 25]
\begin{align}
{{\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,
  \tag{36}
\end{align}
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 diagonal matrix. According to the above decompositions of MAC and BC channel matrixes, the beamforming structures of ${\text{UE}}\;i$ and UAV relay can be proposed as
\begin{align}
{{\mathbf{W}}_i} = {{\mathbf{A}}_{{h_i}}}{{\mathbf{\Lambda }}_i},\;\;\;{{\mathbf{W}}_R} = {{\mathbf{A}}_g}{{\mathbf{\Lambda }}_R}{\mathbf{B}}_h^{ - 1},
  \tag{37}
\end{align}
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 ${\text{UE}}\;i$ and UAV relay, respectively. Plugging (34), (36) and (37) into (4), the receive signal at ${\text{UE}}\;i$ can be written as
\begin{align}
{{\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) + {{\mathbf{n}}_{i,z}}.
  \tag{38}
\end{align}

After using the zero forcing operation for the received signal at ${\text{UE}}\;i$, i.e., ${\mathbf{B}}_{{g_i}}^{ - 1}{{\mathbf{y}}_i}$, we have
\begin{align}
{{\mathbf{\tilde 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}{{{\mathbf{\tilde n}}}_R} + {{{\mathbf{\tilde n}}}_i}} \right) + {{\mathbf{\tilde n}}_{i,z}},
  \tag{39}
\end{align}
where ${{\mathbf{\tilde y}}_i} = {\mathbf{B}}_{{g_i}}^{ - 1}{{\mathbf{y}}_i}$, ${{\mathbf{\tilde n}}_R} = {\mathbf{B}}_h^{ - 1}{{\mathbf{n}}_R}$, ${{\mathbf{\tilde n}}_i} = {\mathbf{B}}_{{g_i}}^{ - 1}{{\mathbf{n}}_i}$ and ${{\mathbf{\tilde 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{align}
  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 \nonumber\\
  {\text{ + }}E_i^{INI} - \frac{T}{2}{\text{tr}}\left( {{{\mathbf{\Lambda }}_i}{\mathbf{\Lambda }}_i^H} \right) \hfill ,
  \tag{40}
\end{align}
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{align}
  &\mathop {\max }\limits_{\left( {{{\mathbf{\Lambda }}_1},\;{{\mathbf{\Lambda }}_2},\;{{\mathbf{\Lambda }}_R},\;{\rho _1},\;{\rho _2}} \right)} \;\;\;E_1^R + E_2^R \hfill
  \tag{41a}
  \\
  &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
  \tag{41b}
  \\
 & \;\;\;\;\;\;{\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
    \tag{41c}
  \\
 & \;\;\;\;\;\;{\text{tr}}\left( {{{\mathbf{\Lambda }}_i}{\mathbf{\Lambda }}_i^H} \right) \leqslant {P_{i,\;Max}},i \in \left\{ {1,2} \right\}, \hfill
    \tag{41d}
  \\
 & \;\;\;\;\;\;0 < {\rho _i} < 1\;,i \in \left\{ {1,2} \right\}. \hfill
  \tag{41e}
\end{align}

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{align}
  E_i^R \geqslant \overset{\lower0.5em\hbox{$\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) \hfill \nonumber\\
  {\text{ + }}E_i^{INI} - \frac{T}{2}{\text{tr}}\left( {{{\mathbf{\Lambda }}_i}{\mathbf{\Lambda }}_i^H} \right) \hfill ,
  \tag{42}
\end{align}
where ${{\mathbf{\Lambda }}_{{g_i}}}$ and ${{\mathbf{\Lambda }}_h}$ are two diagonal matrixes that consist of $N$ diagonal elements of ${{\mathbf{\Psi }}_{{g_i}}}$ and ${{\mathbf{\Psi }}_h}$. For arbitrary positive semi-definite matrixes ${\mathbf{A}}$ and ${\mathbf{X}}$, it has ${\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  is a diagonal matrix that comprises the diagonal elements of ${\mathbf{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{align}
  &\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\nonumber \\
  &\;\; \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 ,
  \tag{43}
\end{align}
\begin{align}
{\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) .
  \tag{44}
\end{align}

With (42), (43) and (44), we can reformulate the problem (41) as
\begin{align}
 & \mathop {\max }\limits_{\left( {{{\mathbf{\Lambda }}_1},\;{{\mathbf{\Lambda }}_2},\;{{\mathbf{\Lambda }}_R},\;{\rho _1},\;{\rho _2}} \right)} \;\;\;\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{E} _1^R + \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{E} _2^R \hfill
  \tag{45a}
  \\
  &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
  \tag{45b}
  \\
 & \;\;\;\;\;\;{\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
  \tag{45c}
  \\
 & \;\;\;\;\;\;{\text{tr}}\left( {{{\mathbf{\Lambda }}_i}{\mathbf{\Lambda }}_i^H} \right) \leqslant {P_{i,\;Max}},i \in \left\{ {1,2} \right\}, \hfill
  \tag{45d}
  \\
 & \;\;\;\;\;\;0 < {\rho _i} < 1\;,i \in \left\{ {1,2} \right\} \hfill .
\tag{45e}
\end{align}

\subsection{Optimize UE power allocation}
In the following, we adopt the AO-based optimization algorithm to solve the intractable and nonconvex problem (45). For given ${{\mathbf{\Lambda }}_R}$, ${\rho _1}$ and ${\rho _2}$, the objective function of the sum retained energy can be rewritten as
\begin{align}
\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{E} _1^R + \overset{\lower0.5em\hbox{$\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,
  \tag{46}
\end{align}
where
\begin{align}
{{\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},
  \tag{47}
\end{align}
\begin{align}
{{\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}
  \tag{48}
\end{align}
and
\begin{align}
  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 \nonumber\\
  & + \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) + E_1^{INI} + E_2^{INI} \hfill.
  \tag{49}
\end{align}

Define ${\lambda _{i,n}}$, ${\lambda _{R,n}}$, ${\lambda _{{h_i},n}}$, ${\lambda _{g,n}}$, ${\lambda _{{\Psi _h},n}}$        and ${\lambda _{{\Psi _{{g_i}}},n}}$ as the n-th diagonal entries of ${{\mathbf{\Lambda }}_i}$, ${{\mathbf{\Lambda }}_R}$, ${{\mathbf{\Sigma }}_i}$, ${{\mathbf{\Sigma }}_g}$, ${{\mathbf{\Lambda }}_h}$ and ${{\mathbf{\Lambda }}_{{g_i}}}$, respectively. Following this definitions, the sum retained energy can be transformed in a scalar form as
\begin{align}
\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{E} _1^R + \overset{\lower0.5em\hbox{$\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,
  \tag{50}
\end{align}
where ${\lambda _{{\Omega _i},n}}$ denotes the n-th diagonal entries of ${{\mathbf{\Omega }}_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{align}
  &\mathop {\max }\limits_{\left( {{{\mathbf{\lambda }}_1},\;{{\mathbf{\lambda }}_2}} \right)} \;\;\;\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
  \tag{51a}
  \\
 & 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
  \tag{51b}
  \\
  &\;\;\;\;\;\;\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,
  \tag{51c}
  \\
 & \;\;\;\;\;\;\sum\limits_{n = 1}^N {\lambda _{i,n}^2}  \leqslant {P_{i,\;Max}},i \in \left\{ {1,2} \right\} \hfill ,
  \tag{51d}
\end{align}
where ${{\mathbf{\lambda }}_i} = {\left[ {\lambda _{i,1}^2, \cdots ,\lambda _{i,N}^2} \right]^T}$ indicates the allocation power of the subchannels at ${\text{UE}}\;i$. Since the linearity of the objective function with the terms of ${{\mathbf{\lambda }}_1}$ and ${{\mathbf{\lambda }}_2}$, the optimization problem (51) is convex and can be readily solved through the classical optimization scheme, i.e., CVX tool.

\subsection{Optimize UAV relay power allocation}
Similarly, for given ${{\mathbf{\Lambda }}_1}$, ${{\mathbf{\Lambda }}_2}$, ${\rho _1}$ and ${\rho _2}$, we have the sum retained energy function with UAV relay power allocation matrix ${{\mathbf{\Lambda }}_R}$ as
\begin{align}
\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{E} _1^R + \overset{\lower0.5em\hbox{$\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,
  \tag{52}
\end{align}
where
\begin{align}
{{\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),
  \tag{53}
\end{align}
\begin{align}
{{\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)
  \tag{54}
\end{align}
and
\begin{align}
\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} - \frac{T}{2}{\text{tr}}\left( {{{\mathbf{\Lambda }}_1}{\mathbf{\Lambda }}_1^H{\text{ + }}{{\mathbf{\Lambda }}_2}{\mathbf{\Lambda }}_2^H} \right).
  \tag{55}
\end{align}

According to (50),  (52) can be transformed in a scalar form as
\begin{align}
\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{E} _1^R + \overset{\lower0.5em\hbox{$\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,
  \tag{56}
\end{align}
where ${\lambda _{{\Omega _3},n}}$ and ${\lambda _{{\Omega _4},n}}$ strand for the n-th diagonal entries of ${{\mathbf{\Omega }}_3}$ and ${{\mathbf{\Omega }}_4}$, respectively. Combining (45) with (56), the original optimization problem in (9) can be restated with scalar variables as
\begin{align}
  &\mathop {\max }\limits_{\left( {\;{{\mathbf{\lambda }}_R}} \right)} \;\;\;\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 \tag{57a}
  \\
 & 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
 \nonumber \\
  &\; \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\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
  \tag{57b}
  \\
  &\;\;\;\;\;\;\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 ,
  \tag{57c}
\end{align}
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 ${{\mathbf{\lambda }}_R}$, the problem (57) is also convex and can be efficiently solved by CVX tool.


\subsection{Optimize PS factors}
For given ${{\mathbf{\Lambda }}_1}$, ${{\mathbf{\Lambda }}_2}$ and ${{\mathbf{\Lambda }}_R}$, the objective function of maximizing the sum retained energy with ${\rho _1}$ and ${\rho _2}$ can be written as
\begin{align}
\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{E} _1^R + \overset{\lower0.5em\hbox{$\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,
  \tag{58}
\end{align}
where
\begin{align}
{\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),
  \tag{59}
\end{align}
\begin{align}
{\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)
  \tag{60}
\end{align}
and
\begin{align}
\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).
  \tag{61}
\end{align}

Based on the SNR constraint in (45), we have
\begin{align}
{\rho _i} \geqslant \frac{{{\gamma _i}{\text{tr}}\left( {{\mathbf{\Lambda }}_{{g_i}}^{ - 1}} \right)\sigma _{i,z}^2}}{{{{\hat \mu }_i}}},
  \tag{62}
\end{align}
where
\begin{align}
  {{\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 \nonumber\\
  & - {\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.
  \tag{63}
\end{align}

Combining (58) and (62), the original optimization problem can be reformulated as
\begin{align}
  &\mathop {\max }\limits_{\left( {{\rho _1},\;{\rho _2}} \right)} \;\;\;{{\hat \beta }_1}\left( {1 - {\rho _1}} \right){\text{ + }}{{\hat \beta }_2}\left( {1 - {\rho _2}} \right){\text{ + }}\hat C \hfill
  \tag{64a}
   \\
  &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
  \tag{64b}
  \\
  &\;\;\;\;\;\;\;0 < {\rho _i} < 1\;,i \in \left\{ {1,2} \right\}. \hfill
  \tag{64c}
\end{align}

From the relationship of variates in (64), the optimal closed-form expression of PS factors can be expressed as
\begin{align}
{\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\}.
  \tag{65}
\end{align}

In summary, the proposed Low-complexity AO scheme based on GSVD for maximizing the sum remained energy can be summarized as Algorithm 2.
\begin{algorithm}[htb]
\caption{Low-complexity beamforming design for maximizing the sum retained energy}
\label{alg:example}
\begin{algorithmic}[1] %这个1 表示每一行都显示数字
%\REQUIRE The transmit power vector ${\bf{p}}$, the step-size $r$, constant\\
%\quad \; $\delta $ and $\varepsilon $. \\ %算法的输入参数：Input
%\ENSURE The optimal transmit power vector ${{\bf{p}}^*}$.  \\ %算法的输出：Output
\STATE {\bf{Initialization:}}  Set ${\rho _1} = {\rho _2} = 0.5$ and ${{\mathbf{\Lambda }}_1} = {{\mathbf{\Lambda }}_2} = {{\mathbf{\Lambda }}_R} = {{\mathbf{I}}_N}$.
\STATE Adopt GSVD and SVD techniques to decompose $\left[ {{\mathbf{H}}_1^H,\;{\mathbf{H}}_2^H} \right]$ and ${\left[ {{\mathbf{G}}_1^T,\;{\mathbf{G}}_2^T} \right]^T}$, respectively, and then obtain ${{\mathbf{\Sigma }}_1}$, ${{\mathbf{\Sigma }}_2}$, ${{\mathbf{\Sigma }}_g}$, ${{\mathbf{\Lambda }}_h}$, ${{\mathbf{\Lambda }}_{{g_1}}}$ and  ${{\mathbf{\Lambda }}_{{g_2}}}$.
\STATE {\bf{Repeat}\\
\texttt{a) Update the transmission power allocation of UEs ${{\mathbf{\Lambda }}_1}$ and ${{\mathbf{\Lambda }}_2}$ by solving the convex problem (51);\\
b) Update the transmission power allocation of UAV relay ${{\mathbf{\Lambda }}_R}$ by solving the convex problem (57);\\
c) Update the optimal PS factors solution ${\rho _{1,\;opt}}$ and ${\rho _{1,\;opt}}$.}}
\STATE {\bf{Until convergence}
\STATE Plug the solved ${{\mathbf{\Lambda }}_1}$, ${{\mathbf{\Lambda }}_2}$ and ${{\mathbf{\Lambda }}_R}$ back into (37) to obtain the beamforming matrixes ${{\mathbf{W}}_1}$, ${{\mathbf{W}}_2}$ and ${{\mathbf{W}}_R}$}.
    %return返回值
\end{algorithmic}
\end{algorithm}

\emph{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.

\section{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 Figure 2, the UAV relay hovers at a constant altitude that is set as ${d_v} = 8$ meters and the horizontal distance between ${\text{UE 1}}$ and ${\text{UE 2}}$ is fixed as ${d_h} = 30$ meters. Moreover, we denote the horizontal distances of ${\text{UE 1}}$-to-relay and relay-to-${\text{UE 1}}$ as $\xi {d_h}$ and $\left( {1 - \xi } \right){d_h}$, respectively. $\xi  \in \left( {0,\;1} \right)$ is the location ratio of UAV relay. It is intuitive that the real distances of ${\text{UE 1}}$-to-relay and relay-to-${\text{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], we model the MIMO channel matrixes as ${{\mathbf{H}}_i} = {{\mathbf{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{H} }}_i}\sqrt {L_i^{ - \varepsilon }} $ and ${{\mathbf{G}}_i} = {{\mathbf{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{G} }}_i}\sqrt {L_i^{ - \varepsilon }} $, where $L_i^{ - \varepsilon }$ strands for the large-scale fading, ${{\mathbf{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{H} }}_i}$ and ${{\mathbf{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{G} }}_i}$ indicate the small-scale Rayleigh fading. Similar to [17], the path loss exponent of the high altitude transmission is set as $\varepsilon  = 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 $\eta  = 0.8$, $\xi  = 0.5$, $N = 2$, $T = 1$ second, ${\rho _i} = 0.5$, ${\gamma _i} = 35$ dB, $E_i^{INI} = 50$ mJ, ${P_{i,\;Max}} = 20$ dBm and ${P_{R,\;Max}} = 20$ dBm unless otherwise explained.

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.

Figure 4 depicts the various performance results of different schemes versus the maximum transmission power of UAV relay. It is intuitive from Figure 4 (a) that both the proposed AO scheme and low-complexity GSVD scheme significantly outperform the average power assignment (APA) scheme in [27]. 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.

In Figure 4 (b), 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 ${P_{R,Max}}$. Moreover, as expected, the cost energy of APA scheme decreases as ${P_{R,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 ${P_{R,Max}}$ grows.

In Figure 4 (c), we investigate the impact of ${P_{R,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 ${P_{R,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 negligible due to the similar iteration procedure.

Figure 5 illustrates the influence of the location ratio of UAV relay $\xi $ 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].

\section{Conclusion}
In this paper, the joint transceiver beamforming and PS factors optimization of UAV-enabled two-way relaying system has 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 sun 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 sun remained energy.














\begin{backmatter}

\section*{Abbreviations}
UAVs:Unmanned aerial vehicles; SWIPT:Simultaneous wireless information and power transfer; UEs:User equipments; DF:Decode-and-forward; PS:Power splitting; TS:Time switching; MIMO:Multiple-input multiple-output; GSVD:Generalized singular value decomposition; EH:Energy harvesting; ID:Information decoding; RF:Radio-frequency; NOMA: Non-orthogonal multiple access; CSI: Channel state information; AF:Amplify-and-forward; MAC:Multiple access; BC:Broadcast; SDP:Semidefinite programming; AO:Alternating optimization

\section*{Acknowledgements}
The authors would like to thank the Editor and anonymous reviewers for their helpful comments and suggestions in improving the quality of this paper. This work was supported by the National Natural Science Foundation of China under Grant 61671184 and 61901137 and in part by Department of Product Development Unit 4G5G Ericsson (China) Communications Co.LTD.

\section*{Authors’ contributions}
JW is the main writer of this paper and proposed the main idea. GW conceived and designed the experiments. GC and BL participated in the design of the study and performed the statistical analysis. RZ and RZ undertook revision works of the paper. All authors read and approved the final manuscript.

\section*{Availability of data and materials}
Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

\section*{Funding}
This work was supported by the National Natural Science Foundation of China under Grant 61671184 and 61901137 and in part by Department of Product Development Unit 4G5G Ericsson (China) Communications Co.LTD.

\section*{Competing interests}
  The authors declare that they have no competing interests.





\begin{thebibliography}{99}
\bibitem{1}
Mu, G., Joint beamforming and power allocation for wireless powered UAV-assisted cooperative NOMA systems. EURASIP J. Wirel. Commun. Networking, 2020(1), 1-14, doi: 10.1007/13638.1687-1499.
\bibitem{2}
N. Zhao et al., UAV-Assisted Emergency Networks in Disasters. IEEE Wireless Commun., vol. 26, no. 1, pp. 45-51, Feb. 2019, doi: 10.1109/MWC.2018.1800160.
\bibitem{3}
Yan, Y., Shi, W., Zhang, X. Design of UAV wireless power transmission system based on coupling coil structure optimization. EURASIP J. Wirel. Commun. Networking, 2020(1), 1-13. doi: 10.1186/s13638-020-01679-4.
\bibitem{4}
N. Zhao et al., Joint Trajectory and Precoding Optimization for UAV-Assisted NOMA Networks. IEEE Trans. Commun., vol. 67, no. 5, pp. 3723-3735, May 2019, doi: 10.1109/TCOMM.2019.2895831.
\bibitem{5}
S. Jia, J. Zhang, H. Zhao and R. Zhang, Relay Selection for Improved Security in Cognitive Relay Networks With Jamming. IEEE Wireless Commun. Lett., vol. 6, no. 5, pp. 662-665, Oct. 2017, doi: 10.1109/LWC.2017.2728528.
\bibitem{6}
C. Li, S. Zhang, P. Liu, F. Sun, J. M. Cioffi and L. Yang, Overhearing Protocol Design Exploiting Intercell Interference in Cooperative Green Networks. IEEE Trans. Veh. Technol., vol. 65, no. 1, pp. 441-446, Jan. 2016, doi: 10.1109/TVT.2015.2389826.
\bibitem{7}
Q. Song, F. Zheng, Y. Zeng and J. Zhang, Joint Beamforming and Power Allocation for UAV-Enabled Full-Duplex Relay. IEEE Trans. Veh. Technol., vol. 68, no. 2, pp. 1657-1671, Feb. 2019, doi: 10.1109/TVT.2018.2889349.
\bibitem{8}
B. Ji, Y. Li, D. Cao, C. Li, s. Mumtaz and D. Wang, Secrecy Performance Analysis of UAV Assisted Relay Transmission for Cognitive Network with Energy Harvesting. IEEE Trans. Veh. Technol., doi: 10.1109/TVT.2020.2989297.
\bibitem{9}
X. Hu, K. Wong, K. Yang and Z. Zheng, UAV-Assisted Relaying and Edge Computing: Scheduling and Trajectory Optimization. IEEE Trans. Wirel. Commun., vol. 18, no. 10, pp. 4738-4752, Oct. 2019, doi: 10.1109/TWC.2019.2928539.
\bibitem{10}
L. Li, T. Chang and S. Cai, UAV Positioning and Power Control for Two-Way Wireless Relaying. IEEE Trans. Wirel. Commun., vol. 19, no. 2, pp. 1008-1024, Feb. 2020, doi: 10.1109/TWC.2019.2950301.
\bibitem{11}
B. Li, J. Yang, H. Yang, G. Liu, R. Ma, X. Peng. Decode-and-forward cooperative transmission in wireless sensor networks based on physical-layer network coding. Wireless Networks, 2019, doi: 10.1007/s11276-019-02092-6.
\bibitem{12}
C. Li, P. Liu, C. Zou, F. Sun, J. M. Cioffi and L. Yang, Spectral-Efficient Cellular Communications With Coexistent One- and Two-Hop Transmissions. IEEE Trans. Veh. Technol., vol. 65, no. 8, pp. 6765-6772, Aug. 2016, doi: 10.1109/TVT.2015.2472456.
\bibitem{13}
W. Lu, S. Hu, X. Liu, C. He and Y. Gong, Incentive Mechanism Based Cooperative Spectrum Sharing for OFDM Cognitive IoT Network. IEEE Trans. Netw. Sci. Eng., May 2019 doi: 10.1109/TNSE.2019.2917071.
\bibitem{14}
C. Li, F. Sun, J. M. Cioffi and L. Yang, Energy Efficient MIMO Relay Transmissions via Joint Power Allocations. IEEE Trans. Circuits Syst., II, Exp. Briefs, vol. 61, no. 7, pp. 531-535, July 2014, doi: 10.1109/TCSII.2014.2327317.
\bibitem{15}
Z. Ding, S. M. Perlaza, I. Esnaola and H. V. Poor, Power Allocation Strategies in Energy Harvesting Wireless Cooperative Networks. IEEE Trans. Wirel. Commun., vol. 13, no. 2, pp. 846-860, Feb. 2014, doi: 10.1109/TWC.2013.010213.130484.
\bibitem{16}
Y. Ma, H. Chen, Z. Lin, Y. Li and B. Vucetic, Distributed and Optimal Resource Allocation for Power Beacon-Assisted Wireless-Powered Communications. IEEE Trans. Commun., vol. 63, no. 10, pp. 3569-3583, Oct. 2015, doi: 10.1109/TCOMM.2015.2468215.
\bibitem{17}
Z. Wen, X. Liu, S. Zheng and W. Guo, Joint Source and Relay Design for MIMO Two-Way Relay Networks With SWIPT. IEEE Trans. Veh. Technol., vol. 67, no. 1, pp. 822-826, Jan. 2018, doi: 10.1109/TVT.2017.2727061.
\bibitem{18}
X. Zhou and Q. Li, Energy Efficiency for SWIPT in MIMO Two-Way Amplify-and-Forward Relay Networks. IEEE Trans. Veh. Technol., vol. 67, no. 6, pp. 4910-4924, June 2018, doi: 10.1109/TVT.2018.2819682.
\bibitem{19}
W. Wang, R. Wang, H. Mehrpouyan, N. Zhao and G. Zhang, Beamforming for Simultaneous Wireless Information and Power Transfer in Two-Way Relay Channels. IEEE Access, vol. 5, pp. 9235-9250, 2017, doi: 10.1109/ACCESS.2017.2701830.
\bibitem{20}
R. Zhang and C. K. Ho, MIMO Broadcasting for Simultaneous Wireless Information and Power Transfer. IEEE Trans. Wirel. Commun., vol. 12, no. 5, pp. 1989-2001, May 2013, doi: 10.1109/TWC.2013.031813.120224.
\bibitem{21}
B. Li and Y. Rong, Joint Transceiver Optimization for Wireless Information and Energy Transfer in Nonregenerative MIMO Relay Systems. IEEE Trans. Veh. Technol., vol. 67, no. 9, pp. 8348-8362, Sept. 2018, doi: 10.1109/TVT.2018.2846556.
\bibitem{22}
J. Wang, G. Wang, H. Yang, B. Li and T. Zhang, Joint Transceiver Optimization for Multiuser Multi-Antenna Relay Systems With Energy Harvesting. IEEE Access, vol. 7, pp. 151156-151167, 2019, doi: 10.1109/ACCESS.2019.2947171.
\bibitem{23}
J. Rostampoor, S. M. Razavizadeh and I. Lee, Energy Efficient Precoding Design for SWIPT in MIMO Two-Way Relay Networks. IEEE Trans. Veh. Technol., vol. 66, no. 9, pp. 7888-7896, Sep. 2017, doi: 10.1109/TVT.2017.2681942.
\bibitem{24}
Q. Qi, X. Chen and D. W. K. Ng, Robust Beamforming for NOMA-Based Cellular Massive IoT With SWIPT. IEEE Trans. Signal Process., vol. 68, pp. 211-224, 2020, doi: 10.1109/TSP.2019.2959246.
\bibitem{25}
R. Wang and M. Tao, Joint Source and Relay Precoding Designs for MIMO Two-Way Relaying Based on MSE Criterion. IEEE Trans. Signal Process., vol. 60, no. 3, pp. 1352-1365, Mar. 2012, doi: 10.1109/TSP.2011.2178598.
\bibitem{26}
R. Mo and Y. H. Chew, MMSE-based joint source and relay precoding design for amplify-and-forward MIMO relay networks. IEEE Trans. Wirel. Commun., vol. 8, no. 9, pp. 4668-4676, Sep. 2009, doi: 10.1109/TWC.2009.081294.
\bibitem{27}
G. Chen, J. Wang, G. Wang, Y. Zou and D. Zhao, Residual Energy Optimization for MIMO SWIPT Two-Way Relaying System. in \textit{Proc. IEEE 89th Veh. Technol. Conf. (VTC Spring)}, Honolulu, HI, USA, 2019, pp. 1-5, doi: 10.1109/VTCFall.2019.8891084.
\bibitem{28}
W. Wang et al., Joint Precoding Optimization for Secure SWIPT in UAV-Aided NOMA Networks. IEEE Trans. Commun., doi: 10.1109/TCOMM.2020.2990994.
\bibitem{29}
S. Yin, Y. Zhao, L. Li and F. R. Yu, UAV-Assisted Cooperative Communications With Time-Sharing Information and Power Transfer. IEEE Trans. Veh. Technol., vol. 69, no. 2, pp. 1554-1567, Feb. 2020, doi: 10.1109/TVT.2019.2956167.
\bibitem{30}
X. Sun, W. Yang and Y. Cai, Secure Communication in NOMA-Assisted Millimeter-Wave SWIPT UAV Networks. IEEE Internet Things J., vol. 7, no. 3, pp. 1884-1897, Mar. 2020, doi: 10.1109/JIOT.2019.2957021.


\end{thebibliography}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                  The Bibliography                       %%
%%                                                         %%
%%  Bmc_mathpys.bst  will be used to                       %%
%%  create a .BBL file for submission.                     %%
%%  After submission of the .TEX file,                     %%
%%  you will be prompted to submit your .BBL file.         %%
%%                                                         %%
%%                                                         %%
%%  Note that the displayed Bibliography will not          %%
%%  necessarily be rendered by Latex exactly as specified  %%
%%  in the online Instructions for Authors.                %%
%%                                                         %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% if your bibliography is in bibtex format, use those commands:
%\bibliographystyle{bmc-mathphys} % Style BST file (bmc-mathphys, vancouver, spbasic).
%\bibliography{bmc_article}      % Bibliography file (usually '*.bib' )
% for author-year bibliography (bmc-mathphys or spbasic)
% a) write to bib file (bmc-mathphys only)
% @settings{label, options="nameyear"}
% b) uncomment next line
%\nocite{label}

% or include bibliography directly:
% \begin{thebibliography}
% \bibitem{b1}
% \end{thebibliography}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Figures                       %%
%%                               %%
%% NB: this is for captions and  %%
%% Titles. All graphics must be  %%
%% submitted separately and NOT  %%
%% included in the Tex document  %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
%% Do not use \listoffigures as most will included as separate files





\section*{Figures}
  \begin{figure}[h!]
  \caption{\csentence{}UAV-enabled two-way relay system with PS-based UEs.
      }
      \end{figure}

\begin{figure}[h!]
  \caption{\csentence{}The location of UAV-enable two-way relay system.}
      \end{figure}

\begin{figure}[h!]
  \caption{\csentence{}Sum retained energy versus the number of iterations.}
      \end{figure}

      \begin{figure}[h!]
  \caption{\csentence{}Various performance results of different schemes versus the maximum transmission power of UAV relay (dBm) ${P_{R,Max}}$. (a) Sum remained energy versus ${P_{R,Max}}$; (b) Sum cost energy versus ${P_{R,Max}}$; (c) Sum harvested energy versus ${P_{R,Max}}$.}
      \end{figure}

      \begin{figure}[h!]
  \caption{\csentence{} Sum remained energy versus the location ratio of UAV relay $\xi $.}
      \end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Tables                        %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Use of \listoftables is discouraged.
%%
%\section*{Tables}
%\begin{table}[h!]
%\caption{Sample table title. This is where the description of the table should go.}
%      \begin{tabular}{cccc}
%        \hline
%           & B1  &B2   & B3\\ \hline
%        A1 & 0.1 & 0.2 & 0.3\\
%        A2 & ... & ..  & .\\
%        A3 & ..  & .   & .\\ \hline
%      \end{tabular}
%\end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Additional Files              %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\section*{Additional Files}
%  \subsection*{Additional file 1 --- Sample additional file title}
%    Additional file descriptions text (including details of how to
%    view the file, if it is in a non-standard format or the file extension).  This might
%    refer to a multi-page table or a figure.
%
%  \subsection*{Additional file 2 --- Sample additional file title}
%    Additional file descriptions text.


\end{backmatter}
\end{document}
