Joint Design of Altitude and Channel Statistics Based Energy Beamforming for UAV-Enabled Wireless Energy Transfer

Abstract

In recent years, UAV-enabled wireless energy transfer (WET) has attracted significant attention for its ability to provide ground devices with efficient and stable power by flexibly navigating three-dimensional (3D) space and utilizing favorable line-of-sight (LoS) channels. At the same time, energy beamforming utilizing multiple antennas, in which energy beams are focused toward devices in desirable directions, has been highlighted as a key technology for substantially enhancing radio frequency (RF)-based WET efficiency. Despite its significant utility, energy beamforming has not been studied in the context of UAV-enabled WET system design. In this paper, we propose the joint design of UAV altitude and channel statistics based energy beamforming to minimize the overall charging time required for all energy-harvesting devices (EHDs) to meet their energy demands while reducing the additional resources and costs associated with channel estimation. Unlike previous works, in which only the LoS dominant channel without small-scale fading was considered, we adopt a more general air-to-ground (A2G) Rician fading channel, where the LoS probability as well as the Rician factor is dependent on the UAV altitude. To tackle this highly nonconvex and nonlinear design problem, we first examine the scenario of a single EHD, drawing insights by deriving an optimal energy beamforming solution in closed form. We then devise efficient methods for jointly designing altitude and energy beamforming in scenarios with multiple EHDs. Our numerical results demonstrate that the proposed joint design considerably reduces the overall charging time while significantly lowering the computational complexity compared to conventional methods.

1. Introduction

Radio frequency (RF)-based wireless energy transfer (WET) is a key technology enabling the development of intelligent and self-sustaining Internet of Everything (IoE) networks in the 6G era [1,2,3,4,5,6,7]. It provides the continuous supply of power to wireless devices over the air, such as wearable electronic devices, extended reality devices, and robotics, without the need for frequent battery replacements or wired power lines. Moreover, compared to traditional power systems, RF-based WET can significantly enhance the quality of service for powering devices by adapting to different physical conditions and service requirements while also improving throughput and robustness [1,2,3,4,5,6,7].

In recent years, unmanned aerial vehicles (UAVs) have gained significant attention in various scenarios due to their deployment flexibility, mobility, and cost-effectiveness, leading to their widespread adoption across various applications, including military operations, cargo delivery, disaster management, and communication platforms [8,9,10,11,12,13,14,15]. They provide greater flexibility in system design and operation in wireless networks as the deployment position and path of a UAV can be adjusted, yielding significant advantages such as coverage enhancements [8,9,10,11,12,13,14,15].

Due to its advantages, UAV-enabled wireless energy transfer (WET) has received great attention for providing ground devices with more efficient and stable power compared to conventional WET systems, which use fixed-location energy transmitters [14,15,16,17,18,19,20,21,22,23,24]. In particular, UAVs are able to move flexibly in three-dimensional (3D) space, utilizing favorable line-of-sight (LoS) channels with ground devices. As such, UAVs can serve as a new type of aerial energy transmitter, reducing the transmission distances for powering devices while avoiding obstacles and shadow fading, even in remote areas where conventional fixed-location energy transmitters are not available. Accordingly, UAV-enabled WET is able to overcome energy bottlenecks as well as meet urgent energy demands, extending the operational lifespan of energy-constrained devices, especially in dynamic or hard-to-reach environments such as smart cities, wireless sensor networks, and maritime networks, among others [14,15,16,17,18]. However, effective UAV trajectory design is required to reap potential gains from UAV-enabled WET.

In practical UAV-enabled WET scenarios, an UAV may move closer to a device in order to reduce the transmission distance for improved power efficiency. This may result in the UAV moving farther away from another device, thereby decreasing the overall energy transfer efficiency. To address this challenge, UAV trajectory design has been extensively studied to improve the energy transfer performance for multiple devices [14,15,16,17,18,19,20,21]. Specifically, in a single UAV-enabled WET system, an optimal one-dimensional (1D) trajectory design at a fixed altitude was proposed to maximize the minimum received energy among all devices during a given charging period [16]. For a more general context, the design of a two-dimensional (2D) UAV trajectory at a fixed altitude was studied to optimize the energy transferred to all devices during a given charging period [17], with further investigations in a more practical scenario [18]. Moreover, the design of a 3D UAV trajectory within a specified altitude range was explored to maximize the total received energy across all devices for a given charging period [19]. The joint design of 2D UAV trajectory and orientation of 1D directional antenna array at the UAV with a fixed altitude was also studied [20] and further extended to a structure that includes a 3D directional antenna array at the UAV [21]. Moreover, a multi-UAV-enabled WET system capable of covering a large area was proposed, and effective trajectory designs for multiple UAVs were studied to enhance energy transfer performance across various scenarios [15,22,23,24].

Energy beamforming has been recognized as a promising technology for significantly increasing the energy transfer efficiency of RF-based WET systems, particularly compared to single-antenna omni-directional transmission [2,3,4]. An energy transmitter with multiple antennas can simultaneously focus energy beams toward devices in the desired directions, thereby overcoming high propagation path loss without increasing transmit power or bandwidth. In this regard, various energy beamforming techniques have been proposed for use in terrestrial wireless networks with perfect channel state information (CSI) as well as with imperfect CSI [3,4,5,6,7].

Moreover, due to the advantage of having LoS channels, beamforming techniques combined with resource allocation and optimization have been widely considered to enhance power transfer performance in various UAV-enabled systems [14,25]. Beamforming combined with optimization of placement and resource allocation has been studied to maximize energy efficiency in wireless-powered UAV communication systems with non-orthogonal multiple access (NOMA) [26]. Additionally, hybrid beamforming with resource allocation has been explored in UAV-enabled wireless-powered mobile edge computing networks [27]. Furthermore, the joint optimization of beamforming, transmit power, power-splitting ratio, and UAV trajectory was proposed to enhance communication performance in UAV-enabled relay networks with wireless power transfer [28]. For UAV-enabled wireless-powered communication networks (WPCNs), beamforming techniques have also been explored using a backscattering scheme [29] and reconfigurable intelligent surfaces (RIS) [30], by jointly optimizing time allocation.

Various beamforming techniques for UAV-enabled systems have been proposed in conjunction with resource allocation and optimization. However, existing studies have primarily focused on WPCNs or simultaneous wireless information and power transfer (SWIPT) systems, with an emphasis on optimizing communication performance. To the best of our knowledge, despite its significant utility, an energy beamforming technique that focuses solely on increasing energy transfer efficiency to simultaneously charge multiple devices, rather than on communication performance, has not been fully studied to realize its potential gains in UAV-enabled WET systems.

Motivated by the aforementioned observations, we focus on energy beamforming design to optimize energy transfer efficiency for UAV-enabled WET systems. In practical WET networks, multiple devices have different energy requirements, and their charging times vary based on the amount of energy needed. Therefore, UAV altitude and energy beamforming must be optimized accordingly. To this end, this paper investigates a joint design of UAV altitude and energy beamforming to minimize the overall charging time required for all energy-harvesting devices (EHDs) to meet their energy requirements. Thus, a large number of EHDs are efficiently served simultaneously while avoiding unnecessary energy transfer from the UAV. Our main contributions are summarized as follows:

• We propose the joint design of UAV altitude- and channel statistics-based energy beamforming, where the EHDs’ energy demands are considered in efficiently and simultaneously serving the EHDs while reducing the additional resources and costs associated with obtaining perfect channel state information. In contrast to existing works on UAV-enabled WET, which consider only the LoS dominant channel without small-scale fading, we adopt the more general air-to-ground (A2G) Rician fading channel while also taking into account the characteristic nature of the aerial channel in practical UAV scenarios, where the channel statistics depend on the altitude of the UAV.
• Due to the highly non-convex and non-linear nature of our design problem, we first jointly optimize UAV altitude and energy beamforming in a single-EHD scenario to draw insights. We derive a solution for optimal energy beamforming in closed form, thereby developing an efficient algorithm with low complexity in obtaining the optimal solution.
• We devise an efficient algorithm to jointly optimize the UAV altitude and energy beamforming in a scenario with multiple EHDs by investigating the optimal conditions as well as the dual problem. Motivated by insights from the design for a single-EHD scenario, we also develop an efficient low-complexity method for determining near-optimal altitude and energy beamforming. Moreover, we explore a sub-optimal design by leveraging weighted-sum energy beamforming in closed form with considerably reduced computational complexity.
• The numerical results demonstrate that compared to conventional methods, the proposed methods can significantly reduce the overall charging time while also decreasing the computational complexity.

The rest of this paper is organized as follows. In Section 2, we introduce the system model and formulate our design problem. In Section 3, we jointly optimize UAV altitude and energy beamforming for a single-EHD scenario. In Section 4, we propose the methods for jointly optimizing UAV altitude and energy beamforming in a multiple-EHD scenario. We evaluate our proposed methods in Section 5 and conclude our paper in Section 6.

2. System Model and Problem Formulation

2.1. System Model

As depicted in Figure 1, we consider a UAV-enabled RF-based wireless energy transfer (WET) system where one UAV simultaneously transfers power to U energy-harvesting devices (EHDs) on the ground, such as IoT devices and sensors. The UAV is equipped with N transmit antennas, and the EHDs are equipped with a single antenna. The UAV hovers over the central EHDs with an altitude $h\in [H_{min},H_{max}]$ in order to conserve power for flight, where $H_{min}$ and $H_{max}$ denote the minimum and maximum altitude of the UAV, respectively. In addition, the maximum WET coverage of the UAV on the ground is denoted as $r_{max}$. The horizontal distance between the UAV and u-th EHD, where $u\in \mathcal{U}\triangleq \{1,...,U\}$, is denoted by $r_u\le r_{max}$.

Let ${\mathbf{w}}_b\in {\mathbb{C}}^N$ be the b-th energy beamforming vector transmitted from the UAV to the EHDs where $b\in \mathcal{B}\triangleq \{1,...,B\}$ and B ($\le$N) is the number of energy beams to be determined. Then, the transmitted signal vector of the UAV for the power transfer, denoted by $\mathbf{x}\in {\mathbb{C}}^N$, is given as

$$\begin{array}{c}\hfill \mathbf{x}=\sum _{b=1}^B{\mathbf{w}}_b{s}_b\end{array}$$

(1)

where $s_b$ is the energy-bearing signal with zero mean and unit variance, i.e., $\mathbb{E}\left[\right|s_b{|}^2]=1$, which can have an arbitrary distribution, since $s_b$ does not carry any information [1,2,3]. In this case, the transmit covariance matrix at the UAV, denoted by $\mathbf{X}\in {\mathbb{H}}^N$, becomes

$$\begin{array}{c}\hfill \mathbf{X}=\sum _{b=1}^B{\mathbf{w}}_b{\mathbf{w}}_b^H⪰0\end{array}$$

(2)

where ${\mathbb{H}}^N$ represents the set of N-by-N Hermitian matrices. It is assumed that the UAV has a transmit sum power constraint of P; then, the transmit covariance matrix of (2) must satisfy $Tr\left(\mathbf{X}\right)\le P$.

2.2. Aerial Channel Model

In this study, we adopt an air-to-ground (A2G) Rician fading channel by considering the characteristic nature of the aerial channel in practical scenarios, where the channel statistics depend on the UAV altitude [11,12,13,14,30,31]. A Rician fading channel consists of a deterministic line-of-sight (LoS) component as well as a random multipath component, i.e., Rayleigh-distributed non-LoS (NLoS). Thus, the A2G Rician fading channel model characterizes the aerial channel by taking into account the dominant LoS channel and the multiplicative effect of both large- and the small-scale fading. The channel vector between the UAV and the u-th EHD, denoted by ${\mathbf{h}}_u\in {\mathbb{C}}^{N\times 1}$, is modeled as [14,30,31].

$$\begin{array}{c}\hfill {\mathbf{h}}_u=\sqrt{{\beta }_0{d}_u^{-{\alpha }_u}}\left(\sqrt{\frac{{\mathcal{K}}_u}{{\mathcal{K}}_u+1}}{\mathbf{a}}_u+\sqrt{\frac{1}{{\mathcal{K}}_u+1}}{\mathbf{h}}_u^{\mathrm{NLoS}}\right)\end{array}$$

(3)

where ${\beta }_0$ is the path loss at the reference distance; $d_u$ and ${\alpha }_u$ denote the distance and path loss exponent between the UAV and u-th EHD, respectively. ${\mathbf{a}}_u\in {\mathbb{C}}^{N\times 1}$ is the deterministic LoS component; ${\mathbf{h}}_u^{\mathrm{NLoS}}\in {\mathbb{C}}^{N\times 1}$ is the NLoS component distributed with ${\mathbf{h}}_u^{\mathrm{NLoS}}\sim \mathcal{C}\mathcal{N}(\mathbf{0},{\mathbf{I}}_N)$; ${\mathcal{K}}_u$ denotes the Rician factor which reflects the power ratio of the LoS component and NLoS component. Assuming there is a uniform linear array (ULA) at the UAV, the LoS component is represented by

$$\begin{array}{c}\hfill {\mathbf{a}}_u={\left[1,e^{-j2\pi \frac{d_a}{{\lambda }_c}sin\left({\varphi }_u\right)},\dots ,e^{-j2\pi (N-1)\frac{d_a}{{\lambda }_c}sin\left({\varphi }_u\right)}\right]}^T\end{array}$$

(4)

where ${\lambda }_c$ is the carrier wavelength, $d_a$ is the space between adjacent antennas, and ${\varphi }_u$ is the angle of departure (AoD) [14,30,31].

Since the propagation characteristics of the A2G channel, such as obstacle density and link quality, are affected by the altitude of the UAV [8,9,10,11,12,13,14], we model not only the probability of a link having a LoS component but also the values of the path loss exponent and Rician factor as a function of UAV altitude, as follows. First, the distance and elevation angle (in radians) between the UAV and u-th EHD are respectively given by

$$\begin{array}{cc}\hfill d_u& =\sqrt{r_u^2+h^2},\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill {\theta }_u& =arctan\left({\displaystyle \frac{h}{r_u}}\right).\hfill \end{array}$$

(6)

Then, the LoS probability can be modeled as a function of the elevation angle [11,12,13,14]:

$$\begin{array}{c}\hfill P_{\mathrm{LoS}}\left({\theta }_u\right)=\frac{1}{1+a_{1}exp\left(-b_1\left(\frac{180}{\pi }{\theta }_u-a_1\right)\right)}\end{array}$$

(7)

where $a_1$ and $b_1$ are positive parameters that are determined by environmental characteristics, i.e., suburban, urban, dense urban, and highrise urban [8,9,10]. The value of the path loss exponent is typically proportional to the LoS probability, so ${\alpha }_u$ can be modeled as [11,12,13,14]

$$\begin{array}{c}\hfill {\alpha }_u=a_2{P}_{\mathrm{LoS}}\left({\theta }_u\right)+b_2.\end{array}$$

(8)

In (8), the coefficients of $a_2$ and $b_2$ are given by

$$\begin{array}{cc}\hfill a_2& =\left({\alpha }_{\frac{\pi }{2}}-{\alpha }_0\right)\frac{1+a_{1}exp\left(a_1{b}_1\right)}{a_{1}exp\left(a_1{b}_1\right)}\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill b_2& ={\alpha }_0-\frac{a_2}{1+a_{1}exp\left(a_1{b}_1\right)}\hfill \end{array}$$

(10)

where ${\alpha }_0$ and ${\alpha }_{\frac{\pi }{2}}$ are the path loss exponents of the ground and aerial links, which are determined through measurement [10,11,12]. Here, $a_2<0$ and $b_2>0$, and $a_2\cong {\alpha }_{\frac{\pi }{2}}-{\alpha }_0$ and $b_2\cong {\alpha }_0$, where the approximations come from $P_{\mathrm{LoS}}\left(0\right)\to 0$ and $P_{\mathrm{LoS}}\left(\frac{\pi }{2}\right)\to 1$. In addition, with respect to ${\theta }_u$, ${\alpha }_u$ is a monotonic function decreasing from ${\alpha }_0$ to ${\alpha }_{\frac{\pi }{2}}$.

In addition, the Rician factor depends on the LoS probability and is thus characterized by the elevation angle of the UAV. In this regard, ${\mathcal{K}}_u$ can be modeled in dB units using an exponential function of ${\theta }_u$ [11,13,32], as follows:

$$\begin{array}{c}\hfill {\mathcal{K}}_u^{\mathrm{dB}}=a_{3}exp\left(b_3{\theta }_u\right)\end{array}$$

(11)

where $a_3$ and $b_3$ are environment- and frequency-dependent constant parameters that are given by [11,13,32]

$$\begin{array}{c}\hfill a_3={\mathcal{K}}_0^{\mathrm{dB}}, b_3=\frac{2}{\pi }ln\left(\frac{{\mathcal{K}}_{\frac{\pi }{2}}^{\mathrm{dB}}}{{\mathcal{K}}_0^{\mathrm{dB}}}\right).\end{array}$$

(12)

In (12), ${\mathcal{K}}_0^{\mathrm{dB}}$ and ${\mathcal{K}}_{\frac{\pi }{2}}^{\mathrm{dB}}$ are the Rician factors in dB units when the elevation angles are 0 rad and $\frac{\pi }{2}$ rad, respectively, which are obtained from measurements [32]. In addition, with respect to ${\theta }_u$, ${\mathcal{K}}_u^{\mathrm{dB}}$ is monotonic function increasing from ${\mathcal{K}}_0^{\mathrm{dB}}$ to ${\mathcal{K}}_{\frac{\pi }{2}}^{\mathrm{dB}}$.

2.3. Formulation of the Energy Harvesting Optimization Problem

By denoting ${\tau }_u$ as the charging time of the u-th EHD, the energy harvested at the u-th EHD is given by [1,2,3,4,5]

$$\begin{array}{c}\hfill Q_u={\zeta }_u{\tau }_{u}Tr\left({\mathbf{h}}_u{\mathbf{h}}_u^H\mathbf{X}\right), \forall u\in \mathcal{U}\end{array}$$

(13)

where ${\zeta }_k$ ($\in (0,1]$) is the energy-harvesting efficiency, which is determined by the circuit. Without loss of generality, we set ${\zeta }_1=\cdots ={\zeta }_U=1$ for simplicity. However, any constant efficiency values can be integrated into our energy-harvesting optimization problem.

When the energy requirement from the u-th EHD is represented by $Q_u^{\mathrm{Req}}>0$, the required charging time at the u-th EHD to meet the energy requirements is given by ${\tau }_u={\displaystyle \frac{Q_u^{\mathrm{Req}}}{Tr\left({\mathbf{h}}_u{\mathbf{h}}_u^H\mathbf{X}\right)}}$. On the other hand, obtaining perfect channel state information (CSI) at the UAV for energy beamforming design requires additional resources and costs for channel estimation, which are not expected to be feasible in practice due to the high energy consumption and computational complexity [2,4,5,15]. Therefore, we explore channel statistics only in the joint design of UAV altitude and energy beamforming. In the context of optimization for energy harvesting, we focus on the long-term channel statistics, i.e., $\mathbb{E}\left[{\mathbf{h}}_u{\mathbf{h}}_u^H\right]$, rather than the perfect CSI of ${\mathbf{h}}_u{\mathbf{h}}_u^H$ [4,5,7]. To this end, we consider the expected required charging time at the u-th EHD, which is represented as

$$\begin{array}{c}\hfill {\tau }_u=\frac{Q_u^{\mathrm{Req}}}{Tr\left(\mathbb{E}\left[{\mathbf{h}}_u{\mathbf{h}}_u^H\right]\mathbf{X}\right)}, \forall u\in \mathcal{U}.\end{array}$$

(14)

In consideration of the different energy requirements of EHDs, i.e., ${\left[Q_u^{\mathrm{Req}}\right]}_{u=1}^U$, the overall charging time required for all EHDs to meet their energy demands becomes $max({\tau }_1,...,{\tau }_U)$ [3,6,7]. Our design objective is to jointly optimize UAV altitude and energy beamforming to minimize the overall charging time required to meet the energy demands of all the EHDs by leveraging only channel statistics. Consequently, our design problem is formulated as

$$\begin{array}{cc}\hfill {\mathit{P}}_1:\underset{h,\mathbf{X}}{min}& max({\tau }_1,...,{\tau }_U)\hfill \end{array}$$

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& {\tau }_u=\frac{Q_u^{\mathrm{Req}}}{Tr\left(\mathbb{E}\left[{\mathbf{h}}_u{\mathbf{h}}_u^H\right]\mathbf{X}\right)}, \forall u\in \mathcal{U},\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill & Tr\left(\mathbf{X}\right)\le P,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & \mathbf{X}⪰\mathbf{0},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill & H_{min}\le h\le H_{max}.\hfill \end{array}$$

(18)

In problem ${\mathit{P}}_1$, constraints (16) and (17) follow from the transmit sum power constraint and (2), respectively, and the constraint (18) indicates the UAV altitude range with $h\in [H_{min},H_{max}]$. In practical UAV operations, altitude adjustments should consider regulatory constraints such as maximum allowable altitude, no-fly zones, and airspace restrictions, as well as practical factors like energy consumption and safety management [33,34,35].

Additionally, from (3), $\mathbb{E}\left[{\mathbf{h}}_u{\mathbf{h}}_u^H\right]$ in (15) is expressed as

$$\begin{array}{cc}\hfill \mathbb{E}\left[{\mathbf{h}}_u{\mathbf{h}}_u^H\right]& ={\beta }_0{d}_u^{-{\alpha }_u}\left(\underset{\mathrm{LoS}}{\underbrace{\frac{{\mathcal{K}}_u}{{\mathcal{K}}_u+1}{\mathbf{a}}_u{\mathbf{a}}_u^H}}+\underset{\mathrm{NLoS}}{\underbrace{\frac{1}{{\mathcal{K}}_u+1}{\mathbf{I}}_N}}\right)\hfill \\ \hfill & \triangleq {\mathbf{G}}_u\left(h\right),\hfill \end{array}$$

(19)

where $d_u$, ${\alpha }_u$, and ${\mathcal{K}}_u$ are optimization functions for the variable h as well as the horizontal distance between the UAV and u-th EHD, i.e., $r_u$. In (19), $\frac{{\mathcal{K}}_u}{{\mathcal{K}}_u+1}{\mathbf{a}}_u{\mathbf{a}}_u^H$ and $\frac{1}{{\mathcal{K}}_u+1}{\mathbf{I}}_N$ indicate the LoS and NLoS components in the long-term channel statistics, respectively.

By introducing an auxiliary variable as $t\triangleq max({\tau }_1,...,{\tau }_U)$, the problem ${\mathit{P}}_1$ is equivalently reformulated into the problem ${\mathit{P}}_2$, as follows:

$$\begin{array}{cc}\hfill {\mathit{P}}_2:\underset{h,\mathbf{X},t}{min}& t\hfill \end{array}$$

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& Tr\left({\mathbf{G}}_u\left(h\right)\mathbf{X}\right)\ge \frac{Q_u^{\mathrm{Req}}}{t}, \forall u\in \mathcal{U},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & Tr\left(\mathbf{X}\right)\le P,\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill & \mathbf{X}⪰\mathbf{0},\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & H_{min}\le h\le H_{max}.\hfill \end{array}$$

(23)

${\mathit{P}}_2$ is a non-convex optimization problem, since the constraints of (20) are non-convex with respect to the UAV’s altitude h. More specifically, in constraint (20), ${\mathbf{G}}_u\left(h\right)$ is non-linear as well as very complex with respect to h for all $u\in \mathcal{U}$, because $d_u$, ${\alpha }_u$, and ${\mathcal{K}}_u$ in ${\mathbf{G}}_u\left(h\right)$ are non-linear and complex functions of h. Therefore, obtaining an optimal solution is, generally, extremely difficult. In order to tackle the problem ${\mathit{P}}_2$, we first solve the problem in a single-EHD scenario, to draw insights, and then establish efficient methods for obtaining an optimal solution in a multiple-EHD scenario.

3. Joint Optimization of Altitude and Energy Beamforming for a Single EHD

In this section, we first solve the problem ${\mathit{P}}_2$ in a single-EHD scenario by deriving energy beamforming in closed form. We then devise a method to jointly optimize altitude and energy beamforming. Next, we develop an efficient low-complexity algorithm that yields the near-optimal altitude of the UAV.

3.1. Optimal Design with a Closed-Form Energy Beamforming Solution

For a single-EHD scenario, i.e., $U=1$, we have $t={\tau }_U=\frac{Q_U^{\mathrm{Req}}}{Tr\left({\mathbf{G}}_U\left(h\right)\mathbf{X}\right)}$ in the problem ${\mathit{P}}_2$, so it can be reformulated into the problem ${\mathit{P}}_3$:

$$\begin{array}{cc}\hfill {\mathit{P}}_3:\underset{h,\mathbf{X}}{min}& \frac{Q_U^{\mathrm{Req}}}{Tr\left({\mathbf{G}}_U\left(h\right)\mathbf{X}\right)}\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& Tr\left(\mathbf{X}\right)\le P, \mathbf{X}⪰\mathbf{0},\hfill \\ \hfill & H_{min}\le h\le H_{max}.\hfill \end{array}$$

Since $Tr\left({\mathbf{G}}_U\left(h\right)\mathbf{X}\right)>0$ must be satisfied at an optimal value, the objective function of the problem ${\mathit{P}}_3$ can be recast as

$$\begin{array}{c}\hfill \underset{h,\mathbf{X}}{min} \frac{Q_U^{\mathrm{Req}}}{Tr\left({\mathbf{G}}_U\left(h\right)\mathbf{X}\right)}=\underset{h,\mathbf{X}}{max} Tr\left({\mathbf{G}}_U\left(h\right)\mathbf{X}\right).\end{array}$$

(24)

From (24), the problem ${\mathit{P}}_3$ is equivalent to the problem ${\mathit{P}}_{3^{\prime }}$:

$$\begin{array}{cc}\hfill {\mathit{P}}_{3^{\prime }}:\underset{h,\mathbf{X}}{max}& Tr\left({\mathbf{G}}_U\left(h\right)\mathbf{X}\right)\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& Tr\left(\mathbf{X}\right)\le P, \mathbf{X}⪰\mathbf{0},\hfill \\ \hfill & H_{min}\le h\le H_{max}.\hfill \end{array}$$

It is also difficult to obtain an optimal solution for h and $\mathbf{X}$ by directly solving the problem ${\mathit{P}}_{3^{\prime }}$, because ${\mathbf{G}}_U\left(h\right)$ is a non-convex and complex non-linear function of h. However, an optimal covariance matrix for $\mathbf{X}$ in relation to the problem ${\mathit{P}}_{3^{\prime }}$ is obtained in closed-form solution via the following lemma.

Lemma 1. 

Optimal energy beamforming in the case of a single EHD is

$$\begin{array}{cc}\hfill {\mathbf{X}}^{★}& =P{\tilde{\mathbf{a}}}_U{\tilde{\mathbf{a}}}_U^H\hfill \end{array}$$

(25)

where ${\tilde{\mathbf{a}}}_U={\displaystyle \frac{{\mathbf{a}}_U}{\parallel {\mathbf{a}}_U\parallel }}$. In addition, with a fixed $h\in [H_{min},H_{max}]$, the optimal value of the problem ${\mathit{P}}_{3^{\prime }}$ is given by

$$\begin{array}{cc}\hfill Tr\left({\mathbf{G}}_U\left(h\right){\mathbf{X}}^{★}\right)& =P{\beta }_0{d}_U^{{\displaystyle -{\alpha }_U}}\frac{{\mathcal{K}}_U{\parallel {\mathbf{a}}_U\parallel }^2+1}{{\mathcal{K}}_U+1}\hfill \end{array}$$

(26)

where $d_U$, ${\alpha }_U$, and ${\mathcal{K}}_U$ depend on h. In (26), $\frac{{\mathcal{K}}_U{\parallel {\mathbf{a}}_U\parallel }^2}{{\mathcal{K}}_U+1}$ and $\frac{1}{{\mathcal{K}}_U+1}$ indicate the impact of the LoS and NLoS components, respectively.

Proof of Lemma 1. 

Refer to Appendix A.    □

From Lemma 1, it is verified that optimal energy beamforming for a single-EHD case, i.e., ${\mathbf{X}}^{★}$, depends only on the LoS component, i.e., ${\mathbf{a}}_U$, of the EHD and is not affected by the UAV altitude h. Consequently, by inserting (25) and (26) in Lemma 1 into the problem ${\mathit{P}}_{3^{\prime }}$, the optimal UAV altitude for a single EHD can be obtained by solving the following problem:

$$\begin{array}{cc}\hfill {\mathit{P}}_{3^{\prime }}^h:\underset{h}{max}& P{\beta }_0{d}_U^{-{\alpha }_U}\frac{{\mathcal{K}}_U{\parallel {\mathbf{a}}_U\parallel }^2+1}{{\mathcal{K}}_U+1}\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& H_{min}\le h\le H_{max}.\hfill \end{array}$$

To tackle the problem ${\mathit{P}}_{3^{\prime }}^h$, we omit $P{\beta }_0$, which is a constant term with respect to h, and we define the objective function in the problem ${\mathit{P}}_{3^{\prime }}^h$ as

$$\begin{array}{cc}\hfill f_U\left(h\right)& =d_U^{-{\alpha }_U}\frac{{\mathcal{K}}_U{\parallel {\mathbf{a}}_U\parallel }^2+1}{{\mathcal{K}}_U+1}\hfill \\ \hfill & ={(r_U^2+h^2)}^{-\frac{1}{2}(a_2{P}_{\mathrm{LoS}}\left({\theta }_U\right)+b_2)}\frac{{\mathcal{K}}_U{\parallel {\mathbf{a}}_U\parallel }^2+1}{{\mathcal{K}}_U+1}\hfill \end{array}$$

(27)

where (27) is obtained by substituting (5) and (8) into $d_U$ and ${\alpha }_U$, respectively. In this case, ${\theta }_U$ is the function of the variable h as well as the EHD’s distance $r_U$, while $P_{\mathrm{LoS}}\left({\theta }_U\right)$ and ${\mathcal{K}}_U$ are complex functions with respect to ${\theta }_U$. Therefore, for simplicity of notation, we define the functions ${\mathbf{\Psi }}_U\left(h\right)$ and ${\mathsf{\Omega }}_U\left(h\right)$ such that $P_{\mathrm{LoS}}\left({\theta }_U\right)\triangleq \frac{1}{1+{\mathbf{\Psi }}_U\left(h\right)}$ and ${\mathcal{K}}_U={10}^{0.1{\mathcal{K}}_U^{\mathrm{dB}}}\triangleq {\mathsf{\Omega }}_U\left(h\right)$, as follows:

$$\begin{array}{cc}\hfill {\mathbf{\Psi }}_U\left(h\right)& =a_{1}exp\left(a_1{b}_1\right)exp\left(-{\displaystyle \frac{180}{\pi }}{b}_{1}arctan\left({\displaystyle \frac{h}{r_U}}\right)\right)\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill {\mathsf{\Omega }}_U\left(h\right)& ={10}^{{\displaystyle {\acute{a}}_{3}exp\left(b_{3}arctan\left({\displaystyle \frac{h}{r_U}}\right)\right)}}\hfill \end{array}$$

(29)

where ${\acute{a}}_3=a_3/10$. Then, by inserting (28) and (29) into (27), we have

$$\begin{array}{cc}\hfill f_U\left(h\right)& ={\left(r_U^2+h^2\right)}^{{\displaystyle -\frac{1}{2}\left(\frac{a_2}{{\mathbf{\Psi }}_U\left(h\right)+1}+b_2\right)}}\frac{\parallel {\mathbf{a}}_U{\parallel }^2{\mathsf{\Omega }}_U\left(h\right)+1}{{\mathsf{\Omega }}_U\left(h\right)+1}.\hfill \end{array}$$

(30)

As represented in (30), the objective function of the problem ${\mathit{P}}_{3^{\prime }}^h$ is non-convex, complex, and non-linear with respect to h, making it generally challenging to obtain a closed-form solution. Fortunately, the objective function of $f_U\left(h\right)$ only has one variable h with the range of $h\in [H_{min},H_{max}]$, so we adopt a 1D exhaustive line search method to obtain the optimal h, i.e., $h^{★}$. The detailed procedure to obtain the optimal altitude of $h^{★}$ in the case of a single EHD is described in Algorithm 1. In this case, the expected required charging time becomes $t^{★}=\frac{Q_U^{\mathrm{Req}}}{Tr\left({\mathbf{G}}_U\left(h^{★}\right)\mathbf{X}\right)}$.

Algorithm 1 Algorithm to find an optimal solution for the problem ${\mathit{P}}_{3^{\prime }}^h$

1: Input: Parameters of the channel statistics ($a_1,b_1,a_2,b_2,a_3,c_3$), the minimum and             maximum altitude of $H_{min}$ and $H_{max}$, the sample grid accuracy $\Delta$,             an EHD’s horizontal distance and LoS component, i.e., $r_U$ and ${\mathbf{a}}_U$   2: Set: $H_{min}\le h\le H_{max}$ with the sample grid $\Delta$, the iteration number $Iter=1$   3: for  $h=H_{min},H_{min}+\Delta ,...,H_{max}$  do   4:        Obtain $f_U\left(h\right)$ from (30)   5:        $Iter$←$Iter+1$   6: end for   7: Total number of iterations: $Ite{r}_{tot}$←$Iter$   8: Obtain $h^{★}=arg {\displaystyle \underset{H_{min}\le h\le H_{max}}{max}} f_U\left(h\right)$   9: Output: The optimal altitude of $h^{★}$

3.2. Efficient Near-Optimal Altitude Design with Low Complexity

In this subsection, we further explore the problem ${\mathit{P}}_{3^{\prime }}^h$ to obtain a near-optimal solution with low complexity. The first-order derivative of $f_U\left(h\right)$ is obtained via

$$\begin{array}{c}\hfill \frac{\partial f_U\left(h\right)}{\partial h}=\frac{{\left(r_U^2+h^2\right)}^{-\frac{a_2}{2({\mathbf{\Psi }}_U\left(h\right)+1)}-\frac{b_2}{2}-1}}{{\mathsf{\Omega }}_U\left(h\right)+1}{g}_U\left(h\right)\end{array}$$

(31)

where $g_U\left(h\right)$ is given by

$$\begin{array}{cc}{g}_U\left(h\right)\triangleq \hfill & (\parallel {\mathbf{a}}_U{\parallel }^2-1)ln\left(10\right){\acute{a}}_3{b}_3\frac{r_{U}exp\left(b_{3}arctan\left({\displaystyle \frac{h}{r_U}}\right)\right){\mathsf{\Omega }}_U\left(h\right)}{{\mathsf{\Omega }}_U\left(h\right)+1}\hfill \\ \hfill & -\frac{a_2}{2}\frac{180}{\pi }{b}_1\frac{r_{U}ln(r_U^2+h^2)\left(\parallel {\mathbf{a}}_U{\parallel }^2{\mathsf{\Omega }}_U\left(h\right)+1\right){\mathbf{\Psi }}_U\left(h\right)}{{\left({\mathbf{\Psi }}_U\left(h\right)+1\right)}^2}\hfill \\ \hfill & -\left(\frac{a_2}{{\mathbf{\Psi }}_U\left(h\right)+1}+b_2\right)\left(\parallel {\mathbf{a}}_U{\parallel }^2{\mathsf{\Omega }}_U\left(h\right)+1\right)h.\hfill \end{array}$$

(32)

As shown in (31), $\frac{\partial f_U\left(h\right)}{\partial h}$ is also a very complex non-linear function with respect to $h\in [H_{min},H_{max}]$ to be optimally solved with a closed-form solution. Nevertheless, we can obtain the necessary conditions for optimal solution of the problem ${\mathit{P}}_{3^{\prime }}^h$ as follows. The critical point of the problem ${\mathit{P}}_{3^{\prime }}^h$ must satisfy $\frac{\partial f_U\left(h\right)}{\partial h}=0$. Additionally, $\frac{\partial f_U\left(h\right)}{\partial h}=0$ is equivalent to $g_U\left(h\right)=0$, because it is readily verified that the first term in (31) always satisfies $\frac{{\left(r_U^2+h^2\right)}^{-\frac{a_2}{2({\mathbf{\Psi }}_U\left(h\right)+1)}-\frac{b_2}{2}-1}}{{\mathsf{\Omega }}_U\left(h\right)+1}>0, \forall h\in [H_{min},H_{max}]$. Moreover, the feasible set of the problem ${\mathit{P}}_{3^{\prime }}^h$ is given as a closed range of $H_{min}\le h\le H_{max}$. As a result, the optimal solution to the problem ${\mathit{P}}_{3^{\prime }}^h$ must satisfy either

$$\begin{array}{cc}\hfill & \left(i\right) g_U\left(h^{★}\right)=0,\mathrm{or}\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill & \left(ii\right) h^{★}=H_{min},\mathrm{or}\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill & \left(iii\right) h^{★}=H_{max}.\hfill \end{array}$$

(35)

From (i)–(iii) in (33)–(35), we adopt the golden-section (GS)-based line search method to determine the near-optimal solution of $h^{★}$, which is known to have less computational complexity than the 1D exhaustive line search method. The detailed procedure of the proposed method to obtain a near-optimal solution for the problem ${\mathit{P}}_{3^{\prime }}^h$ with low complexity is described in Algorithm 2 in the next page. Our proposed method is verified in the numerical results, i.e., Section 5.

Algorithm 2 Proposed method to obtain a near-optimal solution for the problem ${\mathit{P}}_{3^{\prime }}^h$

1: Input: Parameters of the channel statistics ($a_1,b_1,a_2,b_2,a_3,c_3$), the minimum and             maximum altitude of $H_{min}$ and $H_{max}$, the error tolerance ${ϵ}_{\mathrm{GS}}$, $\delta =\frac{(\sqrt{5}-1)}{2}$,             an EHD’s horizontal distance and LoS component, i.e., $r_U$ and ${\mathbf{a}}_U$   2: Set: $h_a=H_{min}$ and $h_b=H_{max}$, $h_{s,a}=h_a+(1-\delta )(h_b-h_a)$ and         $h_{s,b}=h_b-(1-\delta )(h_b-h_a)$, the iteration number $Iter=1$   3: Obtain $g_{s,a}=g_U\left(h_{s,a}\right)$ and $g_{s,b}=g_U\left(h_{s,b}\right)$ from (32), respectively   4: while $|h_b-h_a|<{ϵ}_{\mathrm{GS}}$ do   5:     if $g_{s,a}>g_{s,b},$ then   6:           $h_b$←$h_{s,b}$, $h_{s,b}$← $h_{s,a}$, and $g_{s,b}$←$g_{s,a}$   7:           $h_{s,a}$←$h_a+(1-\delta )(h_b-h_a)$   8:           $g_{s,a}$←$g_U\left(h_{s,a}\right)$ from (32)   9:     else 10:           $h_a$←$h_{s,a}$,$h_{s,a}$← $h_{s,b}$, and $g_{s,a}$←$g_{s,b}$ 11:           $h_{s,a}$←$h_b-(1-\delta )(h_b-h_a)$ 12:           $g_{s,b}$←$g_U\left(h_{s,b}\right)$ from (32) 13:     end if 14:     $Iter$←$Iter+1$ 15: end while 16: Total number of iterations: $Ite{r}_{tot}$←$Iter$ 17: Obtain $h_{\mathrm{Prop}}^{★}=(h_a+h_b)/2$ 18: Output: The proposed near-optimal altitude of $h_{\mathrm{Prop}}^{★}$

4. Joint Optimization of Altitude and Energy Beamforming for Multiple EHDs

In this section, we devise an efficient algorithm in the joint design of UAV altitude and energy beamforming by solving the problem ${\mathit{P}}_2$ in a multiple-EHD scenario (i.e., $U\ge 2$) based on the dual problem. Next, we propose an efficient low-complexity algorithm that yields near-optimal altitude and energy beamforming at the UAV. Finally, we further explore a suboptimal design by developing weighted-sum energy beamforming in closed form, where the computational complexity is considerably reduced by not using an SDP solver.

4.1. Optimal Design with Efficient Algorithm

As mentioned in Section 2, it is generally very difficult to determine the joint optimal solution by directly solving the original problem ${\mathit{P}}_2$. However, the problem ${\mathit{P}}_2$ with a given $h\in [H_{min},H_{max}]$ is a convex optimization problem, so an optimal covariance matrix of $\mathbf{X}$ can be obtained by solving the following problem:

$$\begin{array}{cc}\hfill {\mathit{P}}_2^{\mathbf{X}}:\underset{\mathbf{X},t}{min}& t\hfill \end{array}$$

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& Tr\left({\mathbf{G}}_u\left(h\right)\mathbf{X}\right)\ge \frac{Q_u^{\mathrm{Req}}}{t}, \forall u\in \mathcal{U},\hfill \end{array}$$

(36)

$$Tr\left(\mathbf{X}\right)\le P,$$

(37)

$$\mathbf{X}⪰\mathbf{0}.$$

The problem ${\mathit{P}}_2^{\mathbf{X}}$ is a semidefinite programming (SDP) problem and can therefore be solved using well-known SDP solvers such as CVX [3,36,37]. Therefore, in order to obtain the optimal energy beamforming of ${\mathbf{X}}^{★}$ as well as the optimal UAV altitude $h^{★}$ of the problem ${\mathit{P}}_2$, we can adopt a 1D exhaustive line search method for $h\in [H_{min},H_{max}]$ and iteratively solve the problem ${\mathit{P}}_2^{\mathbf{X}}$ using a fixed h.

On the other hand, the SDP problem is commonly solved based on the interior-point method, so the computational complexity mainly depends on the number and size of variables as well as the number of constraints [36,37]. By defining the prescribed accuracy of an SDP solver as ${ϵ}_{\mathrm{SDP}}>0$, the computational complexity to solve ${\mathit{P}}_2^{\mathbf{X}}$ becomes [37]

$$\begin{array}{c}\hfill \mathcal{O}\left(ln(1/{ϵ}_{\mathrm{SDP}})\sqrt{N+U+1}{N}^2\left(N^4+N^{2}U+U\right)\right).\end{array}$$

(38)

In obtaining ${\mathbf{X}}^{★}$ for a given h, the computational complexity increases exponentially, following an order of approximately $6\frac{1}{2}$ as the number of antennas N increases, which results in a burden on the UAV. To overcome this challenge, we devised a low-complexity method to obtain ${\mathbf{X}}^{★}$ by solving the dual problem instead of the primal problem ${\mathit{P}}_2^{\mathbf{X}}$, as follows.

${\mathit{P}}_2^{\mathbf{X}}$ is a convex optimization problem that satisfies Slater’s condition [3,36]. Hence, we consider the Lagrangian function of ${\mathit{P}}_2^{\mathbf{X}}$. By denoting ${\nu }_u\ge 0$ ($\forall u\in \mathcal{U}$) and ${\nu }_P\ge 0$ as the dual variables associated with the constraints of (36) and (37), respectively, the Lagrangian function is then obtained as [3]

$$\begin{array}{cc}\hfill L\left({\left\{{\nu }_u\right\}}_{u=1}^U,{\nu }_P\right)& =\underset{\mathbf{X}⪰\mathbf{0},t}{inf}\left[t+\sum _{u=1}^U{\nu }_u\left(\frac{Q_u^{\mathrm{Req}}}{t}-Tr\left({\mathbf{G}}_u\left(h\right)\mathbf{X}\right)\right)+{\nu }_P\left(Tr\left(\mathbf{X}\right)-P\right)\right]\hfill \\ \hfill & =\underset{\mathbf{X}⪰\mathbf{0}}{inf} Tr\left[\left({\nu }_P\mathbf{I}-\sum _{u=1}^U{\nu }_u{\mathbf{G}}_u\left(h\right)\right)\mathbf{X}\right]+\underset{t}{inf}\left(t+\frac{{\sum }_{u=1}^U{\nu }_u{Q}_u^{\mathrm{Req}}}{t}\right)-{\nu }_{P}P\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill & =2\sqrt{\sum _{u=1}^U{\nu }_u{Q}_u^{\mathrm{Req}}}-{\nu }_{P}P\hfill \end{array}$$

(40)

where (40) follows from the fact that $Tr\left[\left({\nu }_P\mathbf{I}-{\sum }_{u=1}^U{\nu }_u{\mathbf{G}}_u\left(h\right)\right)\mathbf{X}\right]=0$ as well as ${\nu }_P\mathbf{I}-{\sum }_{u=1}^U{\nu }_u{\mathbf{G}}_u\left(h\right)⪰\mathbf{0}$ must be satisfied in (39) to ensure that $\mathbf{X}⪰\mathbf{0}$ is bounded below, and $\partial L\left({\left\{{\nu }_u\right\}}_{u=1}^U,{\nu }_P\right)/\partial t=0$ gives $t=\sqrt{{\sum }_{u=1}^U{\nu }_u{Q}_u^{\mathrm{Req}}}$ in (39). We omit the details, since the derivations are well described in [3].

Consequently, the dual problem associated with ${\mathit{P}}_2^{\mathbf{X}}$ is formulated as

$$\begin{array}{cc}\hfill \mathit{D}.{\mathit{P}}_2^{\mathbf{X}}:\underset{{\left\{{\nu }_u\right\}}_{u=1}^U, {\nu }_P}{maximize}& 2\sqrt{\sum _{u=1}^U{\nu }_u{Q}_u^{\mathrm{Req}}}-{\nu }_{P}P\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& {\nu }_P\mathbf{I}-\sum _{u=1}^U{\nu }_u{\mathbf{G}}_u\left(h\right)⪰\mathbf{0},\hfill \\ \hfill & {\nu }_u\ge 0, \forall u\in \mathcal{U}, {\nu }_P\ge 0.\hfill \end{array}$$

(41)

The dual problem $\mathit{D}.{\mathit{P}}_2^{\mathbf{X}}$ can also be solved using well-known SDP solvers. By denoting ${\left\{{\nu }_u^{★}\right\}}_{u=1}^U$ and ${\nu }_P^{★}$ as the optimal dual variables attained from the dual problem $\mathit{D}.{\mathit{P}}_2^{\mathbf{X}}$, the optimal solution $t^{★}$ of the primal problem ${\mathit{P}}_2^{\mathbf{X}}$ becomes

$$\begin{array}{c}\hfill t^{★}=\sqrt{\sum _{u=1}^U{\nu }_u^{★}{Q}_u^{\mathrm{Req}}}.\end{array}$$

(42)

Here, we note that the computational complexity to solve the dual problem $\mathit{D}.{\mathit{P}}_2^{\mathbf{X}}$ approximates to

$$\begin{array}{c}\hfill \mathcal{O}\left(ln(1/{ϵ}_{\mathrm{SDP}})\sqrt{N+U+1}U\left(N^3+N^{2}U+U^2\right)\right).\end{array}$$

(43)

To compare (38) and (43) in terms of their computational complexity, we define the ratio between them as

$$\begin{array}{cc}\hfill r_{\mathcal{O}}& =\frac{ln(1/{ϵ}_{\mathrm{SDP}})\sqrt{N+U+1}U\left(N^3+N^{2}U+U^2\right)}{ln(1/{ϵ}_{\mathrm{SDP}})\sqrt{N+U+1}{N}^2\left(N^4+N^{2}U+U\right)}\hfill \\ \hfill & =\frac{U\left(N^3+N^{2}U+U^2\right)}{N^2\left(N^4+N^{2}U+U\right)}.\hfill \end{array}$$

(44)

Substituting $U=N(N-1)$ into (44) with some manipulations, we obtain

$$\begin{array}{cc}\hfill r_{\mathcal{O}}{|}_{U=N(N-1)}& =\frac{(N-1)(2N^2-2N+1)}{N^3+(N-1)(N^2+1)}\hfill \\ \hfill & ⪅\frac{(N-1)(2N^2-2N+1)}{N^3-1+(N-1)(N^2+1)}\hfill \\ \hfill & =\frac{2N^2-2N+1}{2N^2+N+2}\hfill \\ \hfill & <1.\hfill \end{array}$$

(45)

From (44) and (45), it is verified that the computational complexity to solve the dual problem $\mathit{D}.{\mathit{P}}_2^{\mathbf{X}}$ becomes much lower than that of the primal problem ${\mathit{P}}_2^{\mathbf{X}}$ when $U\le N(N-1)$, where the number of EHDs, i.e., U, is much smaller than $N(N-1)$ when considering a practical scenario for UAV-enabled WET [15,16,17,18,19,20,21,22,23,24]. Moreover, the gap in performance increases as U decreases from $N(N-1)$. For example, when $U=N$, we obtain $r_{\mathcal{O}}{|}_{U=N}<\frac{2}{N^2}$ from (44). Hence, the optimal altitude $h^{★}$ is determined using the proposed method, in which the dual problem $\mathit{D}.{\mathit{P}}_2^{\mathbf{X}}$ is solved based on applying the 1D exhaustive line search method within the range of $h\in [H_{min},H_{max}]$.

On the other hand, when we solve the dual problem $\mathit{D}.{\mathit{P}}_2^{\mathbf{X}}$ with a fixed h, rather than the primal problem ${\mathit{P}}_2^{\mathbf{X}}$, we cannot obtain an exact solution for ${\mathbf{X}}^{★}$ and only the optimal condition for ${\mathbf{X}}^{★}$ [3], i.e., $Tr\left[\left({\nu }_P^{★}\mathbf{I}-{\sum }_{u=1}^U{\nu }_u^{★}{\mathbf{G}}_u\left(h\right)\right)\mathbf{X}\right]=0$, which is further equivalent to $\left({\nu }_P^{★}\mathbf{I}-{\sum }_{u=1}^U{\nu }_u^{★}{\mathbf{G}}_u\left(h\right)\right)\mathbf{X}=\mathbf{0}$. Although this optimal condition provides the structure for an optimal covariance matrix ${\mathbf{X}}^{★}$ with its rank profile, obtaining a closed-form solution for ${\mathbf{X}}^{★}$ is generally not available and known to be NP-hard [3,38]. Therefore, as we determine the optimal value of $h^{★}$ from the dual problem $\mathit{D}.{\mathit{P}}_2^{\mathbf{X}}$, the optimal covariance matrix ${\mathbf{X}}^{★}$ is then obtained by solving the primal problem ${\mathit{P}}_2^{\mathbf{X}}$ using the given optimal $h^{★}$ value. The detailed procedure of the proposed method to obtain the optimal altitude ${\mathbf{h}}^{★}$ and optimal energy beamforming ${\mathbf{X}}^{★}$ is described in Algorithm 3.

Algorithm 3 Proposed method of obtaining an optimal solution of the problem ${\mathit{P}}_2$

1: Input: Parameters of the channel statistics ($a_1,b_1,a_2,b_2,a_3,c_3$), the minimum and             maximum altitude of $H_{min}$ and $H_{max}$, the sample grid accuracy $\Delta$,             the sum power constraint P, all EHDs’ horizontal distance, LoS component, and             energy requirements, i.e., ${\left[r_u,{\mathbf{a}}_u,Q_u^{\mathrm{Req}}\right]}_{u=1}^U$   2: Set: $H_{min}\le h\le H_{max}$ with the sample grid $\Delta$, the iteration number $Iter=1$   3: for  $h=H_{min},H_{min}+\Delta ,...,H_{max}$  do   4:       Compute ${\left[{\mathbf{G}}_u\left(h\right)\right]}_{u=1}^U$ from (A5)   5:       Obtain $t^{★}\left(h\right)$ from (42) by solving the dual problem $\mathit{D}.{\mathit{P}}_2^{\mathbf{X}}$   6:       $Iter$←$Iter+1$   7: end for   8: Total number of iterations: $Ite{r}_{tot}$←$Iter$   9: Obtain $h^{★}=arg {\displaystyle \underset{H_{min}\le h\le H_{max}}{min}} t^{★}\left(h\right)$ 10: Obtain $X^{★}$ and $t^{★}$ by solving the primal problem ${\mathit{P}}_2^{\mathbf{X}}$ with ${\left[{\mathbf{G}}_u\left(h^{★}\right)\right]}_{u=1}^U$ 11: Output: The optimal altitude $h^{★}$ and energy beamforming ${\mathbf{X}}^{★}$

4.2. Efficient Near-Optimal Altitude and Energy Beamforming Design with Low Complexity

In this subsection, we propose a low-complexity method of obtaining a near-optimal solution by avoiding the 1D exhaustive line search method, which results in high complexity when using an SDP solver. Motivated by insights from the method of obtaining the near-optimal solution to the single-EHD case in Section 3, we devise the golden-section (GS) line search method with search space reduction, which determines the solution much faster. The details are as follows.

At the u-th EHD, from Lemma 1, the optimal energy beamforming to minimize the expected required charging time, denoted by ${\mathbf{X}}_u^{★}$, is given by

$$\begin{array}{cc}\hfill {\mathbf{X}}_u^{★}& =P{\tilde{\mathbf{a}}}_u{\tilde{\mathbf{a}}}_u^H, \forall u\in \mathcal{U}.\hfill \end{array}$$

(46)

Moreover, the near-optimal UAV altitude to minimize the required charging time, denoted by $h_{\mathrm{Prop},u}^{★}$, is obtained from Algorithm 1.

When a UAV adopts an altitude and energy beamforming as the optimal altitude and energy beamforming, respectively, for the u-th EHD, i.e., $h=h_{\mathrm{Prop},u}^{★}$ and $\mathbf{X}={\mathbf{X}}_u^{★}$ at the UAV, the expected required charging time at the k-th EHD ($k\in \mathcal{U}$), denoted by ${\tau }_k^{★,u}(h_{\mathrm{Prop},u}^{★},{\mathbf{X}}_u^{★})$, then becomes

$$\begin{array}{c}\hfill {\tau }_k^{★,u}(h_{\mathrm{Prop},u}^{★},{\mathbf{X}}_u^{★})=\frac{Q_k^{\mathrm{Req}}}{Tr\left({\mathbf{G}}_k\left(h_{\mathrm{Prop},u}^{★}\right){\mathbf{X}}_u^{★}\right)}, \forall k\in \mathcal{U}.\end{array}$$

(47)

In this case, the overall charging time required to meet the energy demands of all EHDs (i.e., $\forall k\in \mathcal{U}$), denoted by $t_{★,u}^{max}$, is expressed as

$$\begin{array}{c}\hfill t_{★,u}^{max}=max\left({\tau }_1^{★,u}(h_{\mathrm{Prop},u}^{★},{\mathbf{X}}_u^{★}),...,{\tau }_U^{★,u}(h_{\mathrm{Prop},u}^{★},{\mathbf{X}}_u^{★})\right), \forall u\in \mathcal{U},\end{array}$$

(48)

which signifies the overall charging time corresponding to the pair $(h_{\mathrm{Prop},u}^{★},{\mathbf{X}}_u^{★})$ at the UAV.

Among all the overall charging times obtained from (48), i.e., $t_{★,1}^{max}$, $t_{★,2}^{max}$,…, and $t_{★,U}^{max}$, we can identify the lowest time among those achieved by all pairs, i.e., $(h_{\mathrm{Prop},1}^{★},{\mathbf{X}}_1^{★})$, $(h_{\mathrm{Prop},2}^{★},{\mathbf{X}}_2^{★})$,…, and $(h_{\mathrm{Prop},U}^{★},{\mathbf{X}}_U^{★})$, as follows:

$$\begin{array}{c}\hfill u_{min,1}=arg \underset{u\in \mathcal{U}}{min}\left(t_{★,1}^{max},...,t_{★,u}^{max},...,t_{★,U}^{max}\right),\end{array}$$

(49)

with its corresponding pair $(h_{\mathrm{Prop},u_{min,1}}^{★},{\mathbf{X}}_{u_{min,1}}^{★})$. From (49), we infer that the overall charging time is most effectively reduced for $h_{\mathrm{Prop},u_{min,1}}^{★}$ compared with the other pairs ${\left[h_{\mathrm{Prop},k}^{★}\right]}_{k\ne u_{min,1}}^U$. Moreover, the second lowest time is also obtained using

$$\begin{array}{c}\hfill u_{min,2}=arg \underset{u\in \mathcal{U}\setminus u_{min,1}}{min}\left(t_{★,1}^{max},...,t_{★,u}^{max},...,t_{★,U}^{max}\right),\end{array}$$

(50)

with its corresponding pair $(h_{\mathrm{Prop},u_{min,2}}^{★},{\mathbf{X}}_{u_{min,2}}^{★})$. From (49) and (50), we determine the effective altitude range between $h_{\mathrm{Prop},u_{min,1}}^{★}$ and $h_{\mathrm{Prop},u_{min,2}}^{★}$, which is represented by

$$\begin{array}{c}\hfill \left[min\left(h_{\mathrm{Prop},u_{min,1}}^{★},h_{\mathrm{Prop},u_{min,2}}^{★}\right),max\left(h_{\mathrm{Prop},u_{min,1}}^{★},h_{\mathrm{Prop},u_{min,2}}^{★}\right)\right].\end{array}$$

(51)

The overall charging time cannot be effectively reduced for altitudes outside this range, which is taken into account in our proposed method, where the search space $[H_{min},H_{max}]$ is reduced.

Next, in determining the optimal altitude, we adopt the GS line search method within $\left[min\left(h_{\mathrm{Prop},u_{min,1}}^{★},h_{\mathrm{Prop},u_{min,2}}^{★}\right),max\left(h_{\mathrm{Prop},u_{min,1}}^{★},h_{\mathrm{Prop},u_{min,2}}^{★}\right)\right]$, rather than the exhaustive 1D line search method, in order to reduce the computational complexity. Similarly to the procedures for Algorithm 3 and Section 4.1, we determine the optimal altitude $h_{\mathrm{Prop}}^{★}$ by solving the dual problem $\mathit{D}.{\mathit{P}}_2^{\mathbf{X}}$. The detailed procedure for the proposed method to determine near-optimal altitude $h_{\mathrm{Prop}}^{★}$ and energy beamforming ${\mathbf{X}}_{\mathrm{Prop}}^{★}$ is described in Algorithm 4.

Algorithm 4 Proposed method to obtain a near-optimal solution for the problem ${\mathit{P}}_2$

1: Input: Parameters of the channel statistics ($a_1,b_1,a_2,b_2,a_3,c_3$), the minimum and             maximum altitude of $H_{min}$ and $H_{max}$, the error tolerance ${ϵ}_{\mathrm{GS}}$, $\delta =\frac{(\sqrt{5}-1)}{2}$,             the sum power constraint P, all EHDs’ horizontal distance, LoS component, and             energy requirements, i.e., ${\left[r_u,{\mathbf{a}}_u,Q_u^{\mathrm{Req}}\right]}_{u=1}^U$   2: for  $u=1:U$ do   3:       Obtain $h_{\mathrm{Prop},u}^{★}$ from the Algorithm 2   4:       Compute ${\mathbf{X}}_u^{★}=P{\tilde{\mathbf{a}}}_u{\tilde{\mathbf{a}}}_u^H$ from (46)   5:       Compute $t_{★,u}^{max}=max\left({\tau }_1^{★,u}(h_{\mathrm{Prop},u}^{★},{\mathbf{X}}_u^{★}),...,{\tau }_U^{★,u}(h_{\mathrm{Prop},u}^{★},{\mathbf{X}}_u^{★})\right)$ from (48)   6: end for   7: Obtain: $u_{min,1}=arg {\displaystyle \underset{u\in \mathcal{U}}{min}} t_{★,u}^{max}$  and then  $u_{min,2}=arg {\displaystyle \underset{u\in \mathcal{U}\setminus u_{min,1}}{min}} t_{★,u}^{max}$, respectively   8: Set: $h_a=min\left(h_{\mathrm{Prop},u_{min,1}}^{★},h_{\mathrm{Prop},u_{min,2}}^{★}\right)$ and $h_b=max\left(h_{\mathrm{Prop},u_{min,1}}^{★},h_{\mathrm{Prop},u_{min,2}}^{★}\right)$,       $h_{s,a}=h_a+(1-\delta )(h_b-h_a)$ and $h_{s,b}=h_b-(1-\delta )(h_b-h_a)$, $Iter=1$   9: Obtain: $t_{s,a}=t^{★}\left(h_{s,a}\right)$ and $t_{s,b}=t^{★}\left(h_{s,b}\right)$, respectively, from (42), by solving             the dual problem $\mathit{D}.{\mathit{P}}_2^{\mathbf{X}}$ with ${\left[{\mathbf{G}}_u\left(h_{s,a}\right)\right]}_{u=1}^U$ and ${\left[{\mathbf{G}}_u\left(h_{s,b}\right)\right]}_{u=1}^U$, respectively 10: while $|h_b-h_a|<{ϵ}_{\mathrm{GS}}$  do 11:       if $t_{s,a}<t_{s,b},$ then 12:            $h_b$←$h_{s,b}$, $h_{s,b}$← $h_{s,a}$, and $g_{s,b}$←$g_{s,a}$ 13:            $h_{s,a}$←$h_a+(1-\delta )(h_b-h_a)$ 14:            $t_{s,a}$←$t^{★}\left(h_{s,a}\right)$ from (42) by solving the $\mathit{D}.{\mathit{P}}_2^{\mathbf{X}}$ with ${\left[{\mathbf{G}}_u\left(h_{s,a}\right)\right]}_{u=1}^U$ 15:       else 16:            $h_a$←$h_{s,a}$, $h_{s,a}$← $h_{s,b}$, and $g_{s,a}$←$g_{s,b}$ 17:            $h_{s,a}$←$h_b-(1-\delta )(h_b-h_a)$ 18:            $t_{s,b}$←$t^{★}\left(h_{s,b}\right)$ from (42) by solving the $\mathit{D}.{\mathit{P}}_2^{\mathbf{X}}$ with ${\left[{\mathbf{G}}_u\left(h_{s,b}\right)\right]}_{u=1}^U$ 19:       end if 20:       $Iter$←$Iter+1$ 21: end while 22: Total number of iterations: $Ite{r}_{tot}$←$Iter$ 23: Obtain $h_{\mathrm{Prop}}^{★}=(h_a+h_b)/2$ 24: Obtain ${\mathbf{X}}_{\mathrm{Prop}}^{★}$ and $t_{\mathrm{Prop}}^{★}$ by solving the primal problem ${\mathit{P}}_2^{\mathbf{X}}$ with ${\left[{\mathbf{G}}_u\left(h_{\mathrm{Prop}}^{★}\right)\right]}_{u=1}^U$ 25: Output: The proposed near-optimal altitude $h_{\mathrm{Prop}}^{★}$ and energy beamforming ${\mathbf{X}}_{\mathrm{Prop}}^{★}$

4.3. Sub-Optimal Design with Weighted-Sum Energy Beamforming in a Closed-Form Solution

In this subsection, motivated by the insights from Lemma 1, we further devise a solution for weighted-sum energy beamforming, which is obtained in closed form, thereby avoiding the computational complexity introduced by the SDP solver in Algorithm 4.

Since an optimal energy beam for a single EHD to minimize its required charging time is given in (46) from Lemma 1, we consider energy beamforming that consists of the weighted-sum of the optimal beam for each EHD, which has the structure ${\mathbf{X}}^{\mathrm{WS}}={\sum }_{u=1}^U{w}_u{\tilde{\mathbf{a}}}_u{\tilde{\mathbf{a}}}_u^H$, where $w_u\ge 0$ ($u\in U$) denotes the energy weight for the u-th EHD. In this case, it is desirable to guarantee a certain fairness between all EHDs to reduce the overall charging time. Hence, we assign a higher energy weight to the EHDs for which longer charging times are anticipated such that all charging times are the same, i.e., ${\tau }_1=\cdots ={\tau }_U$. For a fixed $h\in [H_{min},H_{max}]$, ignoring the energy amount harvested from other beams, the expected required charging time of the u-th EHD obtained from ${\mathbf{X}}^{\mathrm{WS}}$ becomes ${\tau }_u={\displaystyle \frac{Q_u^{\mathrm{Req}}}{w_{u}Tr\left({\mathbf{G}}_u\left(h\right){\tilde{\mathbf{a}}}_u{\tilde{\mathbf{a}}}_u^H\right)}}={\displaystyle \frac{Q_u^{\mathrm{Req}}}{w_u{\beta }_0{f}_u\left(h\right)}}$. Accordingly, by omitting the common constant term ${\beta }_0$, we set the energy weight for the u-th EHD as

$$\begin{array}{c}\hfill w_u=\frac{Q_u^{\mathrm{Req}}}{f_u\left(h\right)}, \forall u\in \mathcal{U}.\end{array}$$

(52)

As a result, for a given UAV altitude h, our proposed expression for weighted-sum energy beamforming that satisfies a transmit sum power constraint, i.e., $Tr\left({\mathbf{X}}_{\mathrm{Prop}}^{\mathrm{WS}}\right)=P$, is

$$\begin{array}{c}\hfill {\mathbf{X}}_{\mathrm{Prop}}^{\mathrm{WS},★}=\frac{P{\sum }_{u=1}^U{w}_u{\tilde{\mathbf{a}}}_u{\tilde{\mathbf{a}}}_u^H}{Tr\left({\sum }_{u=1}^U{w}_u{\tilde{\mathbf{a}}}_u{\tilde{\mathbf{a}}}_u^H\right)}\end{array}$$

(53)

where ${\mathbf{X}}_{\mathrm{Prop}}^{\mathrm{WS}}⪰\mathbf{0}$ is always ensured considering that $w_u>0$. In addition, we determine the optimal UAV altitude for weighted-sum energy beamforming using the golden-section (GS) line search method together with search space reduction, similarly as for Algorithm 4 in Section 4.2. The detailed procedure is described in Algorithm 5.

Algorithm 5 Algorithm for finding a suboptimal solution

1: Lines 1–8 in Algorithm 4   2: Obtain: $t_{s,a}=max\left({\left[{\tau }_u(h_{s,a},{\mathbf{X}}_{\mathrm{Prop}}^{\mathrm{WS}})\right]}_{u=1}^U\right)$ and $t_{s,b}=max\left({\left[{\tau }_u(h_{s,b},{\mathbf{X}}_{\mathrm{Prop}}^{\mathrm{WS}})\right]}_{u=1}^U\right)$,             respectively, from (52) and (53)   3: while  $|h_b-h_a|<{ϵ}_{\mathrm{GS}}$  do   4:       if $t_{s,a}<t_{s,b},$ then   5:            Same as Lines 12–13 in Algorithm 4   6:            $t_{s,a}$←$max\left({\left[{\tau }_u(h_{s,a},{\mathbf{X}}_{\mathrm{Prop}}^{\mathrm{WS}})\right]}_{u=1}^U\right)$ from (52) and (53)   7:       else   8:            Same as Lines 16–17 in Algorithm 4   9:            $t_{s,b}$←$max\left({\left[{\tau }_u(h_{s,b},{\mathbf{X}}_{\mathrm{Prop}}^{\mathrm{WS}})\right]}_{u=1}^U\right)$ from (52) and (53) 10:       end if 11:       $Iter$←$Iter+1$ 12: end while 13: Total number of iterations: $Ite{r}_{tot}$←$Iter$ 14: Obtain $h_{\mathrm{Prop}}^{\mathrm{WS},★}=(h_a+h_b)/2$, and then ${\mathbf{X}}_{\mathrm{Prop}}^{\mathrm{WS},★}$ from (52) and (53) 15: Output: The proposed altitude $h_{\mathrm{Prop}}^{\mathrm{WS},★}$ and weighted-sum energy beamforming ${\mathbf{X}}_{\mathrm{Prop}}^{\mathrm{WS},★}$

5. Numerical Results

In this section, we evaluate the joint design of UAV altitude and energy beamforming along with our proposed methods. The simulation settings were similar to those in [8,9,10,11,12,13,14], and the detailed parameters used to obtain our numerical results are listed in

Table 1. Environment for evaluation.

Description	Symbol	Value
The transmit sum power	P	1 W
The path loss at the reference distance	${\beta }_0$	1
The path loss exponents of ground links	${\alpha }_0$	3.5
The path loss exponents of aerial links	${\alpha }_{\frac{\pi }{2}}$	2
The Rician factor-related parameters	${\mathcal{K}}_0^{\mathrm{dB}}$	5 dB
	${\mathcal{K}}_{\frac{\pi }{2}}^{\mathrm{dB}}$	15 dB
The minimum altitude of a UAV	$H_{min}$	50 m
The maximum altitude of a UAV	$H_{max}$	250 m
The number of antennas of a UAV	N	16

. Presented in Section 5.1 and Section 5.2 are the performance evaluations for single- and multiple-EHD scenarios, respectively.

5.1. Performance Evaluation in a Single-EHD Scenario

In this subsection, we evaluate our proposed method in a single-EHD scenario. For the Rician fading channel, the EHD’s deterministic LoS component ${\mathbf{a}}_U$ is obtained from (4) with $\frac{d_a}{{\lambda }_c}=0.5$ and an AoD of ${\varphi }_U=\frac{\pi }{3}$.

Figure 2 compares the value of the objective function in problem ${\mathit{P}}_{3^{\prime }}^h$ with respect to the UAV’s altitude (i.e., variable h) for different horizontal distances of an EHD, i.e., $r_U=\{50,100,250\}$ m. As shown in Figure 2, the objective function depends on the UAV’s altitude as well as the horizontal distance of an EHD. Specifically, when $r_U=50$ m and $r_U=250$ m, the objective function is monotonically decreasing and increasing, respectively, as the UAV’s altitude increases. On the other hand, it increases and then decreases when $r_U=100$ m.

Figure 3 shows the optimal UAV altitude according to the horizontal distance of the EHD, obtained with the exhaustive 1D line search method (Algorithm 1) and the proposed low-complexity method (Algorithm 2). Moreover,

Table 2. Comparison of the average execution time (seconds).

	Average Execution Time (Seconds)
Exhaustive 1D line searching (Algorithm 1)	1
Proposed method (Algorithm 2)	$7\times {10}^{-5}$

presents the corresponding average execution time to obtain the optimal value, where we set the same value for the sample grid accuracy $\Delta$ and the error tolerance ${ϵ}_{\mathrm{GS}}$, respectively, i.e., $\Delta ={ϵ}_{\mathrm{GS}}={10}^{-5}$, in Algorithms 1 and 2 for performance comparison. As expected, the proposed method yields almost identical results to the 1D exhaustive line search method with a considerably reduced execution time. Specifically, it was verified that the optimal UAV altitude is the minimum altitude for horizontal distances of up to 50 m, and then it increases from 50 m to 250 m for horizontal distances up to 235 m. Beyond that range, it becomes the maximum altitude. This is because the minimum and maximum altitudes of the UAV are 50 m and 250 m, respectively.

5.2. Performance Evaluation in a Multiple-EHD Scenario

In this subsection, we evaluate our proposed methods in a multiple-EHD scenario. For the Rician fading channel, the u-th EHD’s deterministic LoS component ${\mathbf{a}}_u$ is obtained from (4) with $\frac{d_a}{{\lambda }_c}=0.5$. To evaluate the average overall charging time, we take the average of 100 independent realizations, varying the energy requirements ${\left[Q_u^{\mathrm{Req}}\right]}_{u=1}^U$, horizontal distances ${\left[r_u\right]}_{u=1}^U$, and AoD ${\left[{\varphi }_u\right]}_{u=1}^U$ for each realization. For this purpose, the parameters ${\left[Q_u^{\mathrm{Req}}\right]}_{u=1}^U$, ${\left[r_u\right]}_{u=1}^U$, and ${\left[{\varphi }_u\right]}_{u=1}^U$ are uniformly sampled within specified ranges to capture variability. Specifically, the energy requirements ${\left[Q_u^{\mathrm{Req}}\right]}_{u=1}^U$ are chosen from $[1,10]$ mJ, the horizontal distances ${\left[r_u\right]}_{u=1}^U$ are chosen from ${[0,r_{max}]}_{u=1}^U$ where $r_{max}\in [100,150]$ m, and the angles of departure (AoD) ${\left[{\varphi }_u\right]}_{u=1}^U$ are chosen from $[-\pi /2,\pi /2]$ radians. Moreover, our proposed method was evaluated according to the number of EHDs. For performance comparison, we considered the following six methods:

• A weighted sum-power maximization energy beamforming (W-sum EB) with optimal altitude $H_{\mathrm{Opt}}$;
• Optimal energy beamforming (EB) for the UAV altitude of $H_{min}$;
• Optimal EB for the UAV altitude of $H_{max}$;
• Proposed method $\#1$—efficient near-optimal altitude and EB design with low complexity, i.e., Algorithm 4;
• Proposed method $\#2$—suboptimal design with weighted-sum EB in a closed-form solution, i.e., Algorithm 5;
• Optimal design: optimal EB with optimal altitude, i.e., Algorithm 3 (we set the sample grid accuracy to $\Delta =0.5$ m for the 1D exhaustive line search by considering a practical scenario while avoiding high execution times that make it impossible to use).

Notably, we consider the widely adopted W-sum EB for a comparative analysis with other methods in terms of energy beamforming performance. The W-sum EB was proposed to maximize the weighted sum-power transferred to the EHDs by an energy beam directed toward the dominant eigenvector of the composite weighted channel gain matrix of the EHDs [1,2]. The optimal altitude for the W-sum EB is obtained using a 1D exhaustive line search method.

To examine the robustness of our proposed methods across various practical scenarios, we consider both Urban and Dense Urban environments for the A2G Rician fading channel model. In this case, the environment-dependent LoS probability parameters are given by $(a_1, b_1)=(9.61, 0.16)$ for the Urban environment and $(a_1, b_1)=(12.08, 0.11)$ for the Dense Urban environment [8,9]. First, we evaluate the proposed methods in an Urban environment, and then consider a Dense Urban environment.

5.2.1. Performance Evaluation in an Urban Environment

First, we evaluate the various methods when the horizontal distance is 100 m, i.e., $r_{max}=100$ m. Figure 4a,b show the average UAV altitudes and their standard deviations obtained using proposed method $\#1$, proposed method $\#2$, and the optimal design. Using proposed method $\#1$, the near-optimal altitude of the design is obtained considering the entire range. The performance gap between proposed method $\#1$ and the optimal design is less than 4 m. The performance gap of proposed method $\#2$ is larger than that of proposed method $\#1$. In addition, the UAV’s altitude is within $[70,80]$ m, and it increases as the number of EHDs increases. Moreover, Figure 4b demonstrates that the standard deviation is within $[4,15]$ m, and it increases as the number of EHDs increases.

In Figure 5a,b, the average overall charging times and the average execution times of various methods are compared with respect to the number of EHDs. In Figure 5a, the charging time performance for each method is presented, with values of $[279,1062]$ seconds for W-sum EB with $H_{\mathrm{Opt}}$, $[43,173]$ seconds for optimal EB with $H_{min}$, $[79,267]$ seconds for optimal EB with $H_{max}$, $[14.6,53.5]$ seconds for proposed method $\#1$, $[17.6,101.3]$ seconds for proposed method $\#2$, and $[14.2,52.4]$ seconds for the optimal design. As expected, proposed method $\#1$ outperforms the other methods and achieves a near-optimal charging time over the entire range. The performance gap increases as the number of EHDs increases. Specifically, W-sum EB with $H_{\mathrm{Opt}}$ significantly extends the charging time, as the energy beam is not optimized for minimizing charging time. On the other hand, the optimal design yields a large execution time, increasing from 120 s to 180 s as the number of EHDs increases. By contrast, the execution times range from 7 s to 11 s for proposed method $\#1$, a considerable reduction compared to the optimal design. Meanwhile, proposed method $\#2$, which leverages the weighted-sum EB solution in closed form, yields the shortest execution time of below 0.02 s while yielding lower charging times than the other methods with $H_{min}$ and $H_{max}$.

Next, we evaluate the various methods when the horizontal distance is 150 m, i.e., $r_{max}=150$ m. Figure 6a,b show the average UAV altitude and its standard deviation obtained from proposed method $\#1$, proposed method $\#2$, and the optimal design. Proposed method $\#1$ achieves near-optimal performance, with the performance gap staying within 1.5 m. In addition, proposed method $\#1$ performs comparably to proposed method $\#2$, with the performance gap being within 2.5 m. As expected, the UAV altitude increases as the number of EHDs increases. The UAV altitude is within $[108,122]$ m, which is higher than the average altitude when $r_{max}=100$ m, such as in Figure 4a. In addition, Figure 6b exhibits that the standard deviation is within $[8,25]$ m, which is larger than the standard deviation when $r_{max}=100$ m, such as in Figure 4b.

Figure 7a,b show a comparison of the average overall charging time and the average execution time of various methods with respect to the number of EHDs. In Figure 7a, the charging time performance for each method is presented, such as with values of $[627,2461]$ seconds for W-sum EB with $H_{\mathrm{Opt}}$, $[1370,13710]$ seconds for optimal EB with $H_{min}$, $[85,288]$ seconds for optimal EB with $H_{max}$, $[33.1,120]$ seconds for proposed method $\#1$, $[39.8,229]$ seconds for proposed method $\#2$, and $[32.3,119]$ seconds for the optimal design. Proposed method $\#1$ achieves significantly reduced charging time and near-optimal performance, with a considerable gap compared with other methods as the number of EHDs increases. In particular, there is a substantial increase in the charging time with $H_{min}$, which is much lower than the optimal altitude, as seen in Figure 6. Compared to the charging time performance when $r_{max}=100$, as seen in Figure 5a, the overall charging time increases when $r_{max}=150$. For example, the charging time for the optimal design increases from $[14.2,52.4]$ to $[32.3,119]$. Also, it is shown that the conventional EB, i.e., W-sum EB, yields substantial performance degradation compared to the proposed EB. On the other hand, the performance in terms of execution time is similar for $r_{max}=150$ (Figure 7b) and $r_{max}=100$ (Figure 5b). As expected, proposed method $\#1$ yields a much shorter execution time compared to the optimal design, and the shortest execution time for all methods is achieved using proposed method $\#2$.

5.2.2. Performance Evaluation in a Dense Urban Environment

Next, we evaluate our proposed methods in a Dense Urban environment with $r_{max}=100$ m and $r_{max}=150$ m, respectively.

Figure 8a,b exhibit a comparison of the average overall charging time and the average execution time of various methods when $r_{max}=100$ m. It is shown that proposed method $\#1$ achieves a near-optimal charging time over the entire range, while significantly reducing the execution time compared to the optimal design. Moreover, proposed method $\#2$ achieves the shortest execution time among all methods, while outperforming the conventional methods in terms of overall charging time. However, it results in a longer charging time than proposed method $\#1$. In comparison to the charging time performance when $r_{max}=100$ in the Urban environment, as shown in Figure 5a, the overall charging time increases in the Dense Urban environment. The reason is that the LoS probability decreases in the Dense Urban environment compared to the Urban environment. On the other hand, Figure 8b shows that the average execution time in the Dense Urban environment demonstrates similar performance to that in the Urban environment, as depicted in Figure 5b. This result confirms that the proposed methods reduce computational complexity compared to other methods, even in the Dense Urban environment.

We also evaluate the various methods when $r_{max}=150$ m. Proposed method $\#1$ achieves a significantly reduced charging time and yields a much shorter execution time compared to the optimal design. Also, the shortest execution time for all methods is achieved using proposed method $\#2$. However, proposed method $\#2$ is outperformed by the optimal EB with $H_{max}$ in terms of overall charging time when the number of EHDs exceeds 12. Figure 9a also shows that the overall charging time increases in the Dense Urban environment compared to the charging time performance when $r_{max}=150$ in the Urban environment, as depicted in Figure 7a. Notably, Figure 9b shows that the average execution time in the Dense Urban environment demonstrates similar performance to that in the Urban environment, which verifies the robustness of the proposed methods in terms of complexity reduction.

These results verify that our proposed methods significantly outperforms conventional methods in terms of overall charging time, while achieving near-optimal performance and substantially reducing execution time across various practical scenarios.

6. Conclusions

In this paper, we proposed the joint design of a UAV’s altitude- and channel statistics-based energy beamforming in order to minimize the overall charging time required for all EHDs by considering the A2G Rician fading channel. To solve the formulated highly nonconvex and nonlinear problem, we first optimized our design for a single EHD by deriving optimal energy beamforming in closed form, thereby developing a low-complexity algorithm to obtain the optimal altitude. Then, considering multiple EHDs, we developed efficient algorithms for the joint design of altitude and energy beamforming based on the dual problem. We further explored two efficient methods with low complexity, yielding a near-optimal solution driven by insights from the design for a single-EHD scenario, as well as a sub-optimal solution by leveraging closed-form weighted-sum energy beamforming. The numerical results demonstrate that compared to conventional methods, the proposed joint design can be used to substantially reduce both the overall charging time as well as the computational complexity. While the overall charging time increases in the Dense Urban environment compared to the Urban environment due to a lower LoS probability, the average execution time remains similar, highlighting the robustness of the proposed methods in terms of complexity reduction.

Although the proposed methods are based on long-term channel statistics, they can be extended to scenarios with perfect CSI by replacing the long-term channel statistics with instantaneous CSI. In addition, the proposed method can be adapted to different channel models, such as Nakagami fading, by redefining statistical expectations and modifying the algorithm with necessary mathematical derivations to incorporate new fading characteristics. Extending the proposed methods to provide a comprehensive analysis of scalability, especially when dealing with large numbers of EHDs in diverse environmental conditions, remains one of our ongoing research topics. Additionally, the joint design of altitude and energy beamforming, taking into account practical factors such as the UAV’s energy consumption, is also part of our ongoing work for future research. Moreover, future work will focus on enhancing the proposed model to adapt to dynamic environmental changes that may impact the UAV’s ability to maintain optimal altitude, ensuring robustness in real-world scenarios.