Joint Design of Transmit Waveform and Altitude for Unmanned Aerial Vehicle-Enabled Integrated Sensing and Wireless Power Transfer Systems

Abstract

Recently, unmanned aerial vehicle (UAV)-enabled wireless power transfer (WPT) has received great attention as a promising technology for providing stable power to energy-constrained devices by navigating three-dimensional (3D) space, particularly in challenging environments such as maritime networks and smart cities. Additionally, UAV-enabled radar sensing has gained significant attention as a key technology for future 6G networks, as it enables high-accuracy sensing for various applications, such as target detection and tracking, surveillance, and environmental monitoring, as well as autonomous UAV operation. In this regard, we investigated UAV-enabled integrated sensing and wireless power transfer (ISWPT) systems that combine radar sensing and WPT operations on a unified hardware platform, sharing the same spectrum of resources. In order to accurately sense multiple targets and efficiently transfer power to multiple devices at the same time, we propose a method for jointly designing the transmit waveform and UAV altitude, taking into account the fundamental trade-off between radar sensing performance with the desired beam pattern and WPT performance with the desired harvested power of the devices. We first developed an effective method to obtain the optimal waveform and altitude by solving a challenging non-convex optimization problem. Based on this, we developed another efficient, low-complexity method by exploring a novel transmit waveform and optimizing its parameters to reduce computational complexity and thereby lower power consumption in UAVs. The numerical results verify that the proposed method significantly improves both radar sensing and WPT performance, as well as substantially reduces computational complexity.

1. Introduction

Radio frequency (RF)-based wireless power transfer (WPT) is crucial to developing self-sustaining Internet of Everything (IoE) networks in the 6G era, providing continuous power to devices such as wearable, robotic, and extended reality devices without the need for frequent battery replacement or a wired connection [1,2,3,4]. It adapts to various conditions, enhancing transmission capacity and service quality. Additionally, unmanned aerial vehicles (UAVs) have recently gained great attention for their flexibility, mobility, and cost-effectiveness in various fields, including communication, disaster management, and military operations [5,6,7,8,9,10,11]. In this context, UAV-enabled WPT is widely considered a promising technology for efficiently providing power to ground devices, overcoming the limitations of traditional fixed-location power transmitters [12,13,14,15,16,17,18,19,20,21,22]. By leveraging favorable 3D line-of-sight (LoS) channels, UAVs can serve as aerial power transmitters, reducing transmission distance, avoiding obstacles, and effectively accessing remote areas in challenging environments such as maritime networks and smart cities. Therefore, solutions based on UAV-enabled WPT can address urgent power demands or bottlenecks, thereby extending the lifespan of energy-constrained devices.

In practical scenarios, moving the UAV closer to one device may increase the distance from other devices, thereby reducing the overall energy transfer efficiency. Therefore, there has been extensive research on sophisticated UAV trajectory designs capable of enhancing energy transfer efficiency for multiple ground devices [12,13,14,15,16,17,18,19,20,21,22]; for example, studies have focused on one-dimensional (1D) and two-dimensional (2D) trajectories at a fixed altitude for multiple ground devices in single-UAV-enabled WPT systems [12,13], as well as more practical UAV systems [14], and three-dimensional (3D) trajectories within a certain UAV altitude range [15]. In addition, a directional antenna array for single-UAV-enabled WPT systems was introduced to improve the power transfer efficiency based on the joint design of the trajectory and antenna array orientation [16,17]. UAV altitude and energy beamforming were also jointly designed to maximize power transfer efficiency by directing the energy beams toward the devices at the same time [18]. WPT systems, including multiple UAVs, with particular attention paid to the latter’s trajectories, have been studied in various environments in order to enhance WPT coverage while improving power transfer efficiency [19,20,21,22].

Radar sensing has received significant attention as a key technology for future 6G networks since it allows for the provision of high-accuracy sensing services for various applications [23,24,25,26]. Notably, conventional radar sensing and WPT, although promising, have typically been studied and implemented as separate systems. Recently, a new paradigm called integrated sensing and wireless power transfer (ISWPT), which combines radar sensing and WPT operations on a unified hardware platform sharing the same spectrum resources, has been proposed [25], inspired by integrated sensing and communication (ISAC) systems, where radar sensing and communication operations are integrated into one hardware platform. ISWPT systems provide potential advantages in terms of reduced size, cost, and power consumption of the hardware platform, as well as more efficient use of spectrum resources [25,26]. In their initial work, the authors of [25] characterized the fundamental trade-off between radar sensing and WPT operations by optimizing transmit beamforming. Additionally, the authors of [25] studied a transmit beamforming design for near-field ISWPT systems, in which the devices to be charged are located in the near-field region. However, to the best of our knowledge, UAV-enabled ISWPT systems have not yet been studied, despite the significant utility of UAV-enabled radar sensing in various tasks, such as target detection and tracking, surveillance, environmental monitoring, and autonomous UAV operation.

In this study, we investigate UAV-enabled ISWPT systems, where a UAV is equipped with multiple antennas to sense multiple targets while transferring power to multiple devices on the ground at the same time. We propose a method for jointly designing the transmit waveform and UAV altitude in order to enhance both radar sensing performance and wireless power transfer efficiency. The main contributions of our study are summarized as follows:

• We devise a method for jointly designing the transmit waveform and UAV altitude, taking into account the fundamental trade-off between the desired beam pattern for accurate radar sensing and the desired harvested power of the devices for efficient WPT operation.
• As the corresponding design problem results in challenging non-convex optimization, we first explore an optimal solution and then develop a method to find the optimal waveform and altitude by leveraging the Particle Swarm Optimization (PSO) algorithm. Based on the analysis, we further characterize the maximum desired harvested power of the devices for UAV-enabled ISWPT systems.
• In order to lessen power consumption in UAVs by reducing computational complexity, we also develop an efficient, low-complexity method for jointly designing the transmit waveform and UAV altitude. A novel transmit waveform is devised based on insights from the joint optimization problem, and its design parameters, along with the UAV altitude, are jointly optimized by investigating the optimal conditions and developing an efficient algorithm.
• Our numerical results demonstrate that the proposed joint design method considerably enhances both radar sensing and WPT performance while benefiting from substantially decreased computational complexity.

The rest of this paper is organized as follows. In Section 2, we present the system model and introduce our design problem. In Section 3, we describe the steps taken to obtain the optimal design and investigate the maximum desired harvested power of the devices. We propose an efficient method for jointly designing the waveform and UAV altitude in Section 4. We present the evaluation of the proposed method 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 integrated radar sensing and wireless power transfer (RadWPT) system that senses Q radar targets while transferring power to U energy-harvesting devices (EHDs) on the ground. The UAV in the integrated RadWPT system is equipped with N transmitting antennas and hovers at an altitude $h\in [H_{min},H_{max}]$, where $H_{min}$ and $H_{max}$ are the minimum and maximum altitudes, respectively. In addition, the EHDs on the ground are equipped with a single antenna and are located within the coverage area of the UAV at a horizontal distance $r_u$, where $u\in \mathcal{U}\triangleq \{1,\cdots ,U\}$. In the UAV-enabled integrated RadWPT system, the radar signals transmitted by the UAV can be directly utilized to transfer power to the EHDs because the energy-bearing signals for wireless power transfer (WPT) do not need to carry any information [1,2,3,4], in contrast to conventional integrated sensing and communication (ISAC) systems, which require integrated signals for both radar and communication operations. Therefore, the EHDs can harvest energy from the radar probing signals [25,26].

Let $x_n\left(s\right)$ ($s\in 𝒮\triangleq [1,2,\cdots ,S]$) denote the discrete-time baseband waveform transmitted by the UAV’s n-th ($n\in \mathcal{N}\triangleq [1,2,\cdots ,N]$) antenna at time index s, where S is the number of samples of each transmitted waveform [23,24]. Then, the N-dimensional vector that consists of the N transmitted signals of the s-th waveform is given by $\mathbf{x}\left(s\right)={[x_1\left(s\right),x_2\left(s\right),\cdots ,x_N\left(s\right)]}^T$. The covariance matrix of the waveform transmitted by the UAV, which is expressed as $\mathbf{X}\in {\mathbb{H}}^N$, is obtained as follows:

$$\begin{array}{c}\hfill \mathbf{X}=\mathbb{E}\left[\mathbf{x}\left(s\right)\mathbf{x}{\left(s\right)}^H\right]⪰\mathbf{0},\end{array}$$

(1)

where ${\mathbb{H}}^{N\times N}$ is the set of N-by-N Hermitian matrices. Constraining per-antenna power is common practice when determining radar waveforms [23,24,25,26], so the covariance matrix of the transmitted waveform should meet the following constraint:

$$\begin{array}{c}\hfill {\left[\mathbf{X}\right]}_{n,n}={\displaystyle \frac{P_T}{N}},\forall n\in \mathcal{N},\end{array}$$

(2)

where $P_T$ is the total transmit power of the UAV.

The covariance matrix of the transmitted waveform represents the correlation between the signals transmitted from the antenna array in a MIMO radar system. It defines the power distribution of the signals transmitted by each antenna, the orthogonality between the signals, and the spatial distribution. The objective of signal design in MIMO radar is to optimize the covariance matrix to achieve the desired beam pattern or specific performance metrics [23,24].

2.2. MIMO Radar Sensing

For MIMO radar sensing, we assume that the transmitted probing signal is a narrow-band signal and that the propagation path is non-dispersive, i.e., within the line of sight (LoS) [23,24,25,26]. Then, the received baseband signal at an angular direction $\theta$ can be expressed as

$$\begin{array}{c}\hfill y^{\mathrm{Rad}}(s;\theta )={\mathbf{a}}^H\left(\theta \right)\mathbf{x}\left(s\right),\forall s\in 𝒮,\end{array}$$

(3)

where $\mathbf{a}\left(\theta \right)\in {\mathbb{C}}^{N\times 1}$ is the array steering vector of the direction $\theta$. We adopt uniform linear array (ULA) antennas on the UAV, so the array steering vector is given by [11,27]

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

(4)

where ${\lambda }_c$ is the carrier wavelength and $d_a$ is the spacing between adjacent antennas.

From (3), the power of the probing signal at the angular direction $\theta$ is obtained as [23,24]

$$\begin{array}{c}\hfill P^{\mathrm{Rad}}(\theta ;\mathbf{X})={\mathbf{a}}^H\left(\theta \right)\mathbf{X}\mathbf{a}\left(\theta \right).\end{array}$$

(5)

The spatial spectrum in (5), as a function of the generic direction $\theta$, is termed the beam pattern. The primary purpose of MIMO radar is to direct the transmitted beam toward multiple specified directions, which are usually known to the transmitter, thereby allowing for the collection of more information about the targets illuminated by these beams. To this end, transmit beamforming for MIMO radar probing signals aims to optimize the transmit power over sectors of interest by matching the desired transmit beam pattern to focus the spatial power in the target directions [23,24]. To achieve this, the loss function is defined as the weighted sum of the deviations from the desired beam pattern. The details are given below.

Let ${\left[{\theta }_l\right]}_{l=1}^L$ and $\chi \left({\theta }_l\right)$ be the sample angle grids and the desired transmit beam pattern that covers the location sectors of interest, respectively, where L denotes the number of sample angle grids. Then, the MIMO radar waveform is designed to minimize the mean square error (MSE) between the desired and designed beam patterns, which is given as follows [23,24]:

$$\begin{array}{c}\hfill L_{\mathrm{Rad}}(\mathbf{X},\gamma )={\displaystyle \frac{1}{L}}\sum _{l=1}^L{\left(\gamma \chi \left({\theta }_l\right)-P^{\mathrm{Rad}}({\theta }_l;\mathbf{X})\right)}^2,\end{array}$$

(6)

where $\gamma$ is the scaling factor to be optimized. Such a design is called a beam pattern matching design for MIMO radar sensing.

2.3. Wireless Power Transfer with Aerial Channel Model

We denote the channel vector between the UAV and the u-th EHD on the ground as ${\mathbf{h}}_u\in {\mathbb{C}}^{N\times 1}$; for WPT operation, the harvested power of the u-th EHD is given by [1,2,3,4].

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

(7)

where ${\zeta }_u$ ($\in (0,1]$) is the energy-harvesting efficiency of the u-th EHD, which is determined by the RF circuit. Without loss of generality, we assume that ${\zeta }_1=\cdots ={\zeta }_U=1$ for simplicity. We note that any constant value of ${\zeta }_u$ can be incorporated into our design problem.

For the UAV-enabled WPT to the EHDs on the ground, we adopt an air-to-ground (A2G) Rician fading channel to capture the unique characteristics of the UAV aerial channel in practical scenarios [8,9,10,11,18]. This channel type comprises both a deterministic line-of-sight (LoS) component and a Rayleigh-distributed non-line-of-sight (NLoS) component. In addition, given that the propagation characteristics of the A2G channel, such as obstacle density and link quality, vary with the UAV altitude, the probability of having an LoS link, as well as the path loss exponent and the Rician factor, is modeled as a function of the UAV altitude. Consequently, the A2G Rician fading channel model effectively describes the aerial channel by accounting for both the dominant LoS path and the combined effects of large-scale and small-scale fading.

The A2G Rician fading channel vector between the UAV and the u-th EHD is expressed as [11,18,27]

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

(8)

where ${\beta }_0$ denotes the path loss at the reference distance and $d_u$ is the distance between the UAV and the u-th EHD, which is given by

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

(9)

In addition, ${\alpha }_u$ represents the path loss exponent between the UAV and the u-th EHD. The deterministic LoS component is denoted by $\mathbf{a}\left({\varphi }_u\right)$, where ${\varphi }_u$ indicates the u-th EHD’s angle of departure (AoD), and is constructed based on (4) [18,27]. In addition, ${\mathcal{K}}_u$ denotes the Rician factor, which represents the power ratio between the LoS and NLoS components, where the Rayleigh-distributed NLoS component is expressed as follows: ${\mathbf{h}}_u^{\mathrm{NLoS}}\sim \mathcal{C}\mathcal{N}(\mathbf{0},{\mathbf{I}}_N)\in {\mathbb{C}}^{N\times 1}$ [18,27]. Details of each parameter are reported below.

The probability of having an LoS link is modeled as a function of the UAV altitude [5,6,7,8,9,10,11,18] as follows:

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

(10)

where $a_1>0$ and $b_1>0$ are the environment-dependent (i.e., suburban, urban, dense urban, and high-rise urban) parameters. In addition, ${\psi }_u$ denotes the elevation angle (in radians) between the UAV and the u-th EHD, which is expressed as

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

(11)

In addition, the path loss exponent for the u-th EHD is modeled as a function of $P_{\mathrm{LoS}}\left({\psi }_u\right)$ [8,9,10,11,18]; therefore, it is expressed as follows:

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

(12)

where coefficients $a_2$ and $b_2$ are expressed as

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

respectively, where ${\alpha }_0$ and ${\alpha }_{{\displaystyle \frac{\pi }{2}}}$ denote the path loss exponents of the ground links and aerial links, respectively [7,8,9]. Here, we have $a_2<0$ and $b_2>0$, with $a_2\cong {\alpha }_{{\displaystyle \frac{\pi }{2}}}-{\alpha }_0$ and $b_2\cong {\alpha }_0$ due to $P_{\mathrm{LoS}}\left(0\right)\to 0$ and $P_{\mathrm{LoS}}\left({\displaystyle \frac{\pi }{2}}\right)\to 1$.

Moreover, in the A2G Rician fading channel, the Rician factor for the u-th EHD can be modeled as a function of the u-th EHD’s elevation angle, ${\psi }_u$[8,10,18,28]. By denoting $a_3$ and $b_3$ as the environment- and frequency-dependent constant parameters of the Rician factor, respectively, the Rician factor in dB units is expressed as

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

(13)

In (13), $a_3$ and $b_3$ are given as follows:

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

where ${\mathcal{K}}_0^{\mathrm{dB}}$ and ${\mathcal{K}}_{{\displaystyle \frac{\pi }{2}}}^{\mathrm{dB}}$ denote the Rician factors in dB units when the elevation angles are 0 rad and ${\displaystyle \frac{\pi }{2}}$ rad, respectively [28].

Acquiring the instantaneous channel state information of (8) for the UAV is not practical, as it requires additional resources and costs owing to the high computational complexity and energy consumption of channel estimation [2,3,18,19]. To address this issue, in this study, we concentrated on waveform design for WPT by leveraging only the channel statistics of (8), not instantaneous CSI. Based on (7), the expected harvested power of the u-th EHD becomes

$$\begin{array}{cc}\hfill {\overline{P}}_u^{\mathrm{WPT}}& =Tr\left(\mathbb{E}\left[{\mathbf{h}}_u{\mathbf{h}}_u^H\right]\mathbf{X}\right)\hfill \\ \hfill & =Tr\left({\mathbf{G}}_u\left(h\right)\mathbf{X}\right)\forall u\in \mathcal{U}.\hfill \end{array}$$

(14)

where we denote $\mathbb{E}\left[{\mathbf{h}}_u{\mathbf{h}}_u^H\right]$ by ${\mathbf{G}}_u\left(h\right)$, which is expressed as

$$\begin{array}{c}\hfill {\mathbf{G}}_u\left(h\right)={\beta }_0{d}_u^{-{\alpha }_u}\left({\displaystyle \frac{{\mathcal{K}}_u}{{\mathcal{K}}_u+1}}\mathbf{a}\left({\varphi }_u\right){\mathbf{a}}^H\left({\varphi }_u\right)+{\displaystyle \frac{1}{{\mathcal{K}}_u+1}}{\mathbf{I}}_N\right).\end{array}$$

(15)

In (15), we note that the channel statistics of $d_u$, ${\alpha }_u$, and ${\mathcal{K}}_u$ depend on both the UAV altitude (h) and the u-th EHD’s horizontal distance ($r_u$).

Similar to the beam pattern matching design for MIMO radar sensing, the waveform design of the transmit beamforming of WPT signals, as well as UAV altitude, aims to optimize the received power of the EHDs by matching the desired received power. We denote the desired power at the u-th EHD by $P_u^{\mathrm{Des}}>0$; the WPT waveform is designed to minimize the normalized MSE between the desired and expected harvested power associated with all the EHDs, which is defined as follows:

$$\begin{array}{cc}\hfill L_{\mathrm{WPT}}(\mathbf{X},h)& ={\displaystyle \frac{1}{U}}\sum _{u=1}^U{\left({\displaystyle \frac{P_u^{\mathrm{Des}}-Tr\left({\mathbf{G}}_u\left(h\right)\mathbf{X}\right)}{P_u^{\mathrm{Des}}}}\right)}^2\hfill \\ \hfill & ={\displaystyle \frac{1}{U}}\sum _{u=1}^U{\left(1-{\displaystyle \frac{Tr\left({\mathbf{G}}_u\left(h\right)\mathbf{X}\right)}{P_u^{\mathrm{Des}}}}\right)}^2.\hfill \end{array}$$

(16)

2.4. Problem Formulation

Based on (6) and (16), both MIMO radar sensing performance $L_{\mathrm{Rad}}(\mathbf{X},\gamma )$ and WPT performance $L_{\mathrm{WPT}}(\mathbf{X},h)$ depend on the transmit waveform $\mathbf{X}$. Hence, a fundamental trade-off exists between radar sensing and WPT performance in UAV-enabled integrated RadWPT systems, because the optimal $\mathbf{X}$ that minimizes $L_{\mathrm{Rad}}(\mathbf{X},\gamma )$ may not minimize $L_{\mathrm{WPT}}(\mathbf{X},h)$, and vice versa. In addition, $L_{\mathrm{Rad}}(\mathbf{X},\gamma )$ and $L_{\mathrm{WPT}}(\mathbf{X},h)$ are functions of the scaling factor ($\gamma$) and the UAV altitude (h), respectively; therefore, both $\gamma$ and h have an impact on $\mathbf{X}$. Therefore, $\mathbf{X}$, $\gamma$, and h should be optimized by taking into account the fundamental trade-off between radar sensing and WPT performance.

To achieve a flexible trade-off between the two, we devised a weighted optimization to jointly design the transmit waveform ($\mathbf{X}$), the scaling factor ($\gamma$), and the UAV altitude (h) for UAV-enabled integrated RadWPT systems. As a result, our design problem is formulated as follows:

$$\begin{array}{cc}\hfill {\mathbf{P}}_1:\underset{\mathbf{X},\gamma ,h}{min}& \rho L_{\mathrm{Rad}}(\mathbf{X},\gamma )+(1-\rho )L_{\mathrm{WPT}}(\mathbf{X},h)\hfill \end{array}$$

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& {\left[\mathbf{X}\right]}_{n,n}={\displaystyle \frac{P_T}{N}},\forall n\in \mathcal{N}, \hfill \end{array}$$

(17)

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

(18)

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

(19)

where $0\le \rho \le 1$ is a given weighting factor that determines the performance priority between radar sensing and WPT. In problem ${\mathit{P}}_1$, constraints (17) result from the per-antenna power constraint (2), constraint (1) is the positive semidefinite constraint from (17), and constraint (19) indicates the range of UAV altitude. Note that $\rho =0$ and $\rho =1$ correspond to the scenarios where only WPT or only MIMO radar sensing, respectively, is considered.

3. Optimal Design of Transmit Waveform and UAV Altitude

3.1. Optimal Solution to Our Design Problem

Problem ${\mathit{P}}_1$ is non-convex because the objective function is non-convex with respect to h. More specifically, ${\mathbf{G}}_u\left(h\right)$ for all $u\in \mathcal{U}$ in the objective function is a very complex, non-linear function with respect to h. Thus, it is extremely hard to find an optimal solution in general. To tackle this problem, we first mathematically investigated it and then developed an effective method to solve it, as presented in the next subsection.

For a given $h\in [H_{min},H_{max}]$, $L_{\mathrm{WPT}}(\mathbf{X},h)$ in (16) is a convex function, so the objective function of problem ${\mathit{P}}_1$ is also a convex function. Hence, we can find an optimal waveform ($\mathbf{X}$) and an optimal scaling factor ($\gamma$) by solving the following semidefinite programming (SDP) problem:

$$\begin{array}{cc}\hfill {\mathit{P}}_2^{\forall h}:\underset{\mathbf{X},\gamma }{min}& \rho L_{\mathrm{Rad}}(\mathbf{X},\gamma )+(1-\rho )L_{\mathrm{WPT}}(\mathbf{X},h)\hfill \end{array}$$

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& {\left[\mathbf{X}\right]}_{n,n}={\displaystyle \frac{P_T}{N}},\forall n\in \mathcal{N},\hfill \end{array}$$

(20)

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

(21)

We can numerically obtain the optimal solutions to problem ${\mathit{P}}_2^{\forall h}$ with a well-known SDP solver, such as CVX [1,29,30], and denote the solutions by ${\mathbf{X}}^{★}\left(h\right)$ and ${\gamma }^{★}\left(h\right)$ with a given h.

On the other hand, for a given $\mathbf{X}$ and $\gamma$, $L_{\mathrm{Rad}}(\mathbf{X},\gamma )$ in (16) is a constant term; thus, an optimal altitude (h) can be obtained by solving the following problem:

$$\begin{array}{cc}\hfill {\mathit{P}}_3^{\forall \mathbf{X},\gamma }:\underset{h}{min}& L_{\mathrm{WPT}}(\mathbf{X},h)\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& H_{min}\le h\le H_{max}.\hfill \end{array}$$

In this case, the expected harvested power of the u-th EHD can be expressed as

$$\begin{array}{cc}\hfill Tr\left({\mathbf{G}}_u\left(h\right)\mathbf{X}\right)& ={\beta }_0{d}_u^{-{\alpha }_u}\left({\displaystyle \frac{{\mathcal{K}}_u}{{\mathcal{K}}_u+1}}Tr\left(\mathbf{a}\left({\varphi }_u\right){\mathbf{a}}^H\left({\varphi }_u\right)\mathbf{X}\right)+{\displaystyle \frac{1}{{\mathcal{K}}_u+1}}Tr\left({\mathbf{I}}_N\mathbf{X}\right)\right)\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & ={\beta }_0{\left(r_u^2+h^2\right)}^{{\displaystyle -\frac{1}{2}\left(\frac{a_2}{{\Psi }_{r_u}\left(h\right)+1}+b_2\right)}}{\displaystyle \frac{Tr\left(\mathbf{a}\left({\varphi }_u\right){\mathbf{a}}^H\left({\varphi }_u\right)\mathbf{X}\right){\Omega }_{r_u}\left(h\right)+P_T}{{\Omega }_{r_u}\left(h\right)+1}},\hfill \end{array}$$

(23)

where (23) is obtained by substituting (10)–(13) into (22), and ${\Psi }_{r_u}\left(h\right)$ and ${\Omega }_{r_u}\left(h\right)$ are functions with respect to the optimization variable (h) as well as the u-th EHD’s distance ($r_u$), which are defined as follows:

$$\begin{array}{cc}\hfill {\Psi }_{r_u}\left(h\right)& =a_{1}exp\left(-b_1\left[{\displaystyle \frac{180}{\pi }}arctan\left({\displaystyle \frac{h}{r_u}}\right)-a_1\right]\right),\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {\Omega }_{r_u}\left(h\right)& ={10}^{{\displaystyle 0.1a_{3}exp\left(b_{3}arctan\left({\displaystyle \frac{h}{r_u}}\right)\right)}}.\hfill \end{array}$$

(25)

In (24) and (25), ${\Psi }_{r_u}\left(h\right)$ and ${\Omega }_{r_u}\left(h\right)$ result from $P_{\mathrm{LoS}}\left({\psi }_u\right)\triangleq {\displaystyle \frac{1}{1+{\Psi }_{r_u}\left(h\right)}}$ and ${10}^{0.1{\mathcal{K}}_u^{\mathrm{dB}}}\triangleq {\Omega }_{r_u}\left(h\right)$, respectively, for notation simplicity.

As shown in (23), the objective function of problem ${\mathit{P}}_3^{\forall \mathbf{X},\gamma }$ is a non-convex, very complex, non-linear function of h; thus, finding an optimal h by analytically solving it is extremely difficult. Therefore, we devised a Particle Swarm Optimization (PSO)-based method to solve our design problem, as presented in the next subsection.

3.2. Method for Obtaining the Optimal Solution Based on the PSO Algorithm

Particle Swarm Optimization (PSO) is a meta-heuristic global optimization algorithm inspired by swarm intelligence [31,32,33].

It is easy to implement, making it widely applicable to various optimization problems. It is known for its faster convergence compared to other complex global optimization algorithms, and it can be effectively applied to non-linear or complex functions. Additionally, PSO requires relatively few parameters to be tuned and has simple settings, allowing it to be flexibly applied to different problems [34,35,36]. In a PSO algorithm, a swarm of particles interacts within a given space to find the optimal solution. Each particle has a position and velocity within the solution space and adjusts its movement based on both the best solution it has found and the best solution found by the entire swarm. The algorithm iteratively updates each particle’s position and speed, improves the quality of the path with the fitness function, and achieves the optimal solution with the best fitness value. Thanks to its advantages, the PSO algorithm has been widely adopted to optimize trajectory and path planning for UAVs [31,32,33].

More specifically, for a general D-dimensional optimization problem, a swarm of M particles explores the D-dimensional variable space. The m-th particle has position ${\mathit{P}}_m^j=\left[p_{m,1}^{\left(j\right)},p_{m,2}^{\left(j\right)},\cdots ,p_{m,D}^{\left(j\right)}\right]\in {\mathbb{R}}^D$ and velocity ${\mathbf{v}}_m^{\left(j\right)}=\left[v_{m,1}^{\left(j\right)},v_{m,2}^{\left(j\right)},\cdots ,v_{m,D}^{\left(j\right)}\right]\in {\mathbb{R}}^D$ at the j-th iteration $j=1,2,\cdots ,J$, where J represents the total number of iterations [31,32,33]. The PSO algorithm updates the position of the m-th particle as follows:

$$\begin{array}{c}\hfill {\mathbf{p}}_m^{\left(j\right)}={\mathbf{p}}_m^{(j-1)}+{\mathbf{v}}_m^{\left(j\right)},\end{array}$$

(26)

where ${\mathbf{v}}_m^{(j+1)}$ is the rate of the next movement used to obtain its new position and is updated as follows:

$$\begin{array}{c}\hfill v_{m,d}^{\left(j\right)}=\sigma ·v_{m,d}^{(j-1)}+z_1·r_{m,d,1}^{\left(j\right)}·\left(p_{m,d}^{best}-p_{m,d}^{(j-1)}\right)+z_2·r_{m,d,2}^{\left(j\right)}·\left(g_d^{best}-p_{m,d}^{(j-1)}\right),d=1,2,\cdots ,D,\end{array}$$

(27)

where ${\mathbf{p}}_m^{best}=\left[p_{m,1}^{best},p_{m,2}^{best},\cdots ,p_{m,D}^{best}\right]\in {\mathbb{R}}^D$ is the position corresponding to the best solution that has been found by the m-th particle and ${\mathbf{g}}^{best}=\left[g_1^{best},g_2^{best},\cdots ,g_D^{best}\right]\in {\mathbb{D}}^N$ is the position of the best particle in the whole swarm [31,32,33]. In addition, $\sigma \in (0,1)$ denotes the inertia weight controlling the momentum effect for preserving the velocity, and $z_1$ and $z_2$ are the coefficients controlling the ability to learn from the local components and the whole swarm, respectively. $r_{m,d,1}^{(j+1)}$ and $r_{m,d,2}^{(j+1)}$ are random numbers that are uniformly distributed in the range of $[0,1]$.

Next, we devised a method to obtain the optimal solution to problem ${\mathit{P}}_1$, which consists of finding the optimal UAV altitude (h) based on the PSO algorithm for the one-dimensional problem, i.e., $D=1$, detailed as follows: At the j-th iteration, the m-th particle has altitude $p_m^{\left(j\right)}$ and velocity $v_m^{\left(j\right)}$, where $p_m^{\left(j\right)}$ must satisfy $p_m^{\left(j\right)}\in [H_{min},H_{max}]$. For a given altitude $p_m^{\left(j\right)}$, the optimal waveform ${\mathbf{X}}^{★}\left(p_m^{\left(j\right)}\right)$ and scaling factor ${\gamma }^{★}\left(p_m^{\left(j\right)}\right)$ can be obtained by solving problem ${\mathit{P}}_2^{\forall h}$ with $h=p_m^{\left(j\right)}$. In this case, the fitness function, denoted by $f_{\mathrm{Fit}}\left(p_m^{\left(j\right)}\right)$, becomes the objective function of problem ${\mathit{P}}_2^{\forall h}$, which is written as follows:

$$\begin{array}{c}\hfill f_{\mathrm{Fit}}\left(p_m^{\left(j\right)}\right)=\rho L_{\mathrm{Rad}}\left({\mathbf{X}}^{★}\left(p_m^{\left(j\right)}\right),{\gamma }^{★}\left(p_m^{\left(j\right)}\right)\right)+(1-\rho )L_{\mathrm{WPT}}\left({\mathbf{X}}^{★}\left(p_m^{\left(j\right)}\right),p_m^{\left(j\right)}\right).\end{array}$$

(28)

Then, altitude $p_m^{\left(j\right)}$ and velocity $v_m^{\left(j\right)}$ are updated by minimizing $f_{\mathrm{Fit}}\left(p_m^{\left(j\right)}\right)$. Finally, we obtain the optimal altitude, $h^{★}=g^{best}$, based on the PSO algorithm and then find the optimal $\mathbf{X}\left(g^{best}\right)$ and scaling factor $\gamma \left(g^{best}\right)$ by solving problem ${\mathit{P}}_2^{\forall h}$ with $h=g^{best}$. The detailed procedure is described in Algorithm 1.

In Algorithm 1, the overall computational complexity depends primarily on the SDP solver, which finds an optimal waveform by solving problem ${\mathit{P}}_2^{\forall h}$ for each particle at each iteration. In this case, the SDP solver utilizes the interior-point method, with its computational complexity depending on the number and size of both variables and the number of constraints [29,30]. Specifically, calculating $L_{\mathrm{Rad}}(\mathbf{X},\gamma )$ with sample angle grids ${\left[{\theta }_l\right]}_{l=1}^L$ and $L_{\mathrm{WPT}}(\mathbf{X},h)$ with U EHDs in the objective function is equivalent to having additional L and U linear matrix equality constraints in the problem, respectively. Hence, the overall computational complexity of Algorithm 1 is

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

(29)

where ${ϵ}_{\mathrm{SDP}}>0$ is the prescribed accuracy of the SDP solver. Here, we note that the computational complexity of Algorithm 1, used to obtain the optimal solution, depends not only on the number of antennas (N), the number of EHDs (U), and the number of particles and iterations for the PSO algorithm (i.e., M and J) but also on the number of sample angle grids (L). Typically, L is relatively large for sophisticated beam patterns, as in our case. For example, ${\left[{\theta }_l\right]}_{l=1}^L$ is set as uniform samples in the range of $[-{90}^{\circ },{90}^{\circ }]$ with a resolution of $0.1^{\circ }$ [23,24].

Algorithm 1 Proposed algorithm to obtain an optimal solution to problem ${\mathit{P}}_1$

1: Input: A sampled angle grid ${\left[{\theta }_l\right]}_{l=1}^L$ and a desired transmit beam pattern $\chi \left({\theta }_l\right)$,             parameters on the channel statistics ($a_1,b_1,a_2,b_2,a_3,b_3$), minimum and             maximum altitude of $H_{min}$ and $H_{max}$, sum power constraint $P_T$,             all EHDs’ horizontal distance, LoS component, and desired power,             i.e., ${\left[r_u,{\mathbf{a}}_u,P_u^{\mathrm{Des}}>0\right]}_{u=1}^U$, coefficients for PSO algorithm ($\sigma ,z_1,z_2$) 2: for Particles $m=1,2,\cdots ,M$ do 3:       Initialize $p_m^{\left(0\right)}=H_{min}+r_m^{\left(0\right)}(H_{max}-H_{min})$ and $v_m^{\left(0\right)}=0$ 4:       Initialize $p_m^{best}$←$p_m^{\left(0\right)}$ 5:       Obtain ${\mathbf{X}}^{★}\left(p_m^{best}\right)$ and ${\gamma }^{★}\left(p_m^{best}\right)$ by solving problem ${\mathit{P}}_2^{\forall h}$ with a given $h=p_m^{best}$ 6:       Compute the fitness function value $f_{\mathrm{Fit}}\left(p_m^{best}\right)$ from (28) 7: end for 8: Initialize $g^{best}$←$arg\underset{{\left[p_m^{best}\right]}_{m=1}^M}{\min }{f}_{\mathrm{Fit}}\left(p_m^{best}\right)$ 9: for Iterations $j=1,2,\cdots ,J$ do 10:     for Particles $m=1,2,\cdots ,M$ do 11:           Update $v_m^{\left(j\right)}$ and $p_m^{\left(j\right)}$ from (27) and (26), respectively. 12:           if $p_m^{\left(j\right)}<H_{min}$ then 13:              $p_m^{\left(j\right)}$←$H_{min}$ 14:           else if $p_m^{\left(j\right)}>H_{max}$ then 15:              $p_m^{\left(j\right)}$←$H_{max}$ 16:           end if 17:           Obtain ${\mathbf{X}}^{★}\left(p_m^{\left(j\right)}\right)$ and ${\gamma }^{★}\left(p_m^{\left(j\right)}\right)$ by solving problem ${\mathit{P}}_2^{\forall h}$ with a given $h=p_m^{\left(j\right)}$ 18:           Compute the fitness function value $f_{\mathrm{Fit}}\left(p_m^{\left(j\right)}\right)$ from (28) 19:           if $f_{\mathrm{Fit}}\left(p_m^{\left(j\right)}\right)<f_{\mathrm{Fit}}\left(p_m^{best}\right)$ then 20:              $p_m^{best}$←$p_m^{\left(j\right)}$ 21:           end if 22:       end for 23:       Update $g^{best}$←$arg\underset{{\left[p_m^{best}\right]}_{m=1}^M}{\min }{f}_{\mathrm{Fit}}\left(p_m^{best}\right)$ 24: end for 25: Obtain an optimal altitude $h^{★}=g^{best}$ 26: Obtain the optimal waveform ${\mathbf{X}}^{★}\left(g^{best}\right)$ and scaling factor ${\gamma }^{★}\left(g^{best}\right)$ by solving problem ${\mathit{P}}_2^{\forall h}$ with $h=g^{best}$ 27: Output: Optimal solution to problem ${\mathit{P}}_1$, i.e., ${\mathbf{X}}^{★}$, ${\gamma }^{★}$, and $h^{★}$

3.3. Analysis of Maximum Desired Power of an EHD

In this subsection, we focus on the maximum desired power for the u-th EHD, i.e., the maximum value of $P_u^{\mathrm{Des}}$. When considering only WPT operation with a single EHD, i.e., $\rho =0$ and $U=1$, the expected harvested power of the U-th EHD, i.e., (14), depends on the transmit waveform ($\mathbf{X}$) and the UAV altitude (h). Therefore, the maximum harvested power of the U-th EHD can be obtained by solving the following problem:

$$\begin{array}{cc}\hfill {\mathit{P}}_4:\underset{\mathbf{X},h}{max}& Tr\left({\mathbf{G}}_U\left(h\right)\mathbf{X}\right)\hfill \end{array}$$

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& {\left[\mathbf{X}\right]}_{n,n}={\displaystyle \frac{P_T}{N}},\forall n\in \mathcal{N},\hfill \end{array}$$

(30)

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

(31)

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

(32)

To solve problem ${\mathit{P}}_4$ using mathematical analysis, we first relax the per-antenna power constraints (30) and then consider the sum power constraint, i.e., $Tr\left(\mathbf{X}\right)=P$. Here, (30) is a special case of the sum power constraint. Accordingly, problem ${\mathit{P}}_4$ is relaxed as follows:

$$\begin{array}{cc}\hfill {\mathit{P}}_4^{\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)=P,\mathbf{X}⪰\mathbf{0},\hfill \\ \hfill & H_{min}\le h\le H_{max}.\hfill \end{array}$$

As shown in (23), the objective function of problem ${\mathit{P}}_4^{\prime }$ is a non-convex, very complex, non-linear function of h, so finding optimal solutions for $\mathbf{X}$ and h by directly solving the problem is very difficult. However, the optimal waveform ($\mathbf{X}$) of problem ${\mathit{P}}_4^{\prime }$ is obtained from Lemma 1 in [18] as follows:

$$\begin{array}{cc}\hfill {\mathbf{X}}^{★}& =P_T\tilde{\mathbf{a}}\left({\varphi }_U\right){\tilde{\mathbf{a}}}^H\left({\varphi }_U\right),\hfill \end{array}$$

(33)

where $\tilde{\mathbf{a}}\left({\varphi }_U\right)={\displaystyle \frac{\mathbf{a}\left({\varphi }_U\right)}{\parallel \mathbf{a}\left({\varphi }_U\right)\parallel }}$. In other words, it is verified that optimal beamforming for a single EHD is the well-known maximum-ratio transmission (MRT) beamforming. Intuitively, this result is reasonable because MRT beamforming is known to maximize the received power of a single device in MIMO communication systems.

In this case, by (33), we have

$$\begin{array}{cc}\hfill {\left[{\mathbf{X}}^{★}\right]}_{n,n}& ={\left[P_T\tilde{\mathbf{a}}\left({\varphi }_U\right){\tilde{\mathbf{a}}}^H\left({\varphi }_U\right)\right]}_{n,n}\hfill \\ \hfill & ={\displaystyle \frac{P_T}{\parallel \mathbf{a}\left({\varphi }_U\right){\parallel }^2}}{\left[\mathbf{a}\left({\varphi }_U\right){\mathbf{a}}^H\left({\varphi }_U\right)\right]}_{n,n}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{P_T}{N}}{\left[\mathbf{a}\left({\varphi }_U\right){\mathbf{a}}^H\left({\varphi }_U\right)\right]}_{n,n}\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{P_T}{N}},\hfill \end{array}$$

(35)

where (34) and (35) are due to the fact that $\parallel \mathbf{a}\left({\varphi }_U\right){\parallel }^2=N$ and ${\left[\mathbf{a}\left({\varphi }_U\right){\mathbf{a}}^H\left({\varphi }_U\right)\right]}_{n,n}=1$, respectively, which can be readily obtained from (4).

As a result, the optimal waveform of (33) satisfies the per-antenna power constraints (30); therefore, (33) is the optimal waveform of problem ${\mathit{P}}_4$. By substituting (33) into (23), the objective function of problem ${\mathit{P}}_4$ can be rewritten as follows:

$$\begin{array}{cc}\hfill Tr\left({\mathbf{G}}_U\left(h\right){\mathbf{X}}^{★}\right)& ={\beta }_0{\left(r_U^2+h^2\right)}^{{\displaystyle -\frac{1}{2}\left(\frac{a_2}{{\Psi }_{r_U}\left(h\right)+1}+b_2\right)}}{\displaystyle \frac{P_T\left(N{\Omega }_{r_U}\left(h\right)+1\right)}{{\Omega }_{r_U}\left(h\right)+1}}\hfill \\ \hfill & \triangleq g_U\left(h\right).\hfill \end{array}$$

(36)

Thus, the optimal altitude (h) can be obtained by solving the following problem:

$$\begin{array}{cc}\hfill {\mathit{P}}_4^{″}:\underset{h}{max}& g_U\left(h\right)\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& H_{min}\le h\le H_{max}.\hfill \end{array}$$

Problem ${\mathit{P}}_4^{″}$ has only one variable, h, with $h\in [H_{min},H_{max}]$, so we can easily adopt a 1D exhaustive line search method to obtain the optimal altitude $h^{★}$. Finally, the maximum desired power of an EHD is expressed as $g_U\left(h^{★}\right)$; therefore, we need to have $P_u^{\mathrm{Des}}$ less than $g_u\left(h^{★}\right)$ for all $u\in U$ in WPT operation.

4. Proposed Efficient Design for Transmit Waveform and UAV Altitude

In this section, we investigate the proposed efficient low-complexity design to jointly optimize the transmit waveform and UAV altitude by devising a novel transmit waveform and jointly optimizing the design parameters, the scaling factor, and the UAV altitude.

To decrease the computational complexity derived from the SDP solver, based on the insights from MRT beamforming in (33), we propose a novel waveform for UAV-enabled integrated RadWPT systems, whose structure is as follows:

$$\begin{array}{cc}\hfill {\mathbf{X}}_{\mathrm{Prop}}& =\left(1-\sum _{k=1}^U{w}_k\right){\mathbf{X}}_{\mathrm{Rad}}+P_T\sum _{k=1}^U{w}_k\tilde{\mathbf{a}}\left({\varphi }_k\right){\tilde{\mathbf{a}}}^H\left({\varphi }_k\right),\hfill \end{array}$$

(37)

where $w_k\in [0,1]$ ($\forall k\in \mathcal{U}$), satisfying ${\sum }_{k=1}^U{w}_k\le 1$, is the weighting factor to be optimized; it controls the priority of the k-th EHD’s optimal beamforming. ${\mathbf{X}}_{\mathrm{Rad}}$ is the optimal waveform when considering only MIMO radar sensing, i.e., $\rho =1$, and it can be obtained by solving the following problem:

$$\begin{array}{cc}\hfill {\mathit{P}}_5:\underset{{\mathbf{X}}_{\mathrm{Rad}},\gamma }{min}& \rho L_{\mathrm{Rad}}({\mathbf{X}}_{\mathrm{Rad}},\gamma )\hfill \end{array}$$

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& {\left[{\mathbf{X}}_{\mathrm{Rad}}\right]}_{n,n}={\displaystyle \frac{P_T}{N}},\forall n\in \mathcal{N},\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill & {\mathbf{X}}_{\mathrm{Rad}}⪰\mathbf{0}.\hfill \end{array}$$

(39)

Here, it can be readily verified that (37) satisfies the per-antenna power constraints and the positive semidefinite constraint:

$$\begin{array}{cc}\hfill {\left[{\mathbf{X}}_{\mathrm{Prop}}\right]}_{n,n}& ={\displaystyle \frac{P_T}{N}},\forall n\in \mathcal{N},\hfill \\ \hfill {\mathbf{X}}_{\mathrm{Prop}}& ⪰\mathbf{0},\hfill \end{array}$$

(40)

because ${\left[{\mathbf{X}}_{\mathrm{Rad}}\right]}_{n,n}={\displaystyle \frac{P_T}{N}}$ and ${\mathbf{X}}_{\mathrm{Rad}}⪰\mathbf{0}$, ${\left[P_T\tilde{\mathbf{a}}\left({\varphi }_k\right){\tilde{\mathbf{a}}}^H\left({\varphi }_k\right)\right]}_{n,n}={\displaystyle \frac{P_T}{N}}$ and $\tilde{\mathbf{a}}\left({\varphi }_k\right){\tilde{\mathbf{a}}}^H\left({\varphi }_k\right)⪰\mathbf{0}$, and $w_k\in [0,1]$ ($\forall k\in \mathcal{U}$), with ${\sum }_{k=1}^U{w}_k\le 1$. In our waveform design, the weighting factors should be determined to minimize the objective function of problem ${\mathit{P}}_1$.

First, by substituting ${\mathbf{X}}_{\mathrm{Prop}}$ into (5) for MIMO radar sensing operation, we have

$$\begin{array}{cc}\hfill P^{\mathrm{Rad}}({\theta }_l;{\mathbf{X}}_{\mathrm{Prop}})& ={\mathbf{a}}^H\left({\theta }_l\right){\mathbf{X}}_{\mathrm{Prop}}\mathbf{a}\left({\theta }_l\right)\hfill \\ \hfill & ={\mathbf{a}}^H\left({\theta }_l\right){\mathbf{X}}_{\mathrm{Rad}}\mathbf{a}\left({\theta }_l\right)+\sum _{k=1}^U{w}_k{\mathbf{a}}^H\left({\theta }_l\right)\left\{\tilde{\mathbf{a}}\left({\varphi }_k\right){\tilde{\mathbf{a}}}^H\left({\varphi }_k\right)-{\mathbf{X}}_{\mathrm{Rad}}\right\}\mathbf{a}\left({\theta }_l\right)\hfill \\ \hfill & =b_{\mathrm{Rad},l}+{\mathbf{c}}_{\mathrm{Rad},l}^T\mathbf{w},\forall l\in [1,2,\cdots ,L]\triangleq \mathcal{L},\hfill \end{array}$$

(41)

where $\mathbf{w}={[w_1,\cdots w_U]}^T$, and $b_{\mathrm{Rad},l}\in \mathbb{R}$ and ${\mathbf{c}}_{\mathrm{Rad},l}\in {\mathbb{R}}^{U\times 1}$ ($\forall l\in \mathcal{L}$) are expressed as

$$\begin{array}{cc}\hfill b_{\mathrm{Rad},l}& ={\mathbf{a}}^H\left({\theta }_l\right){\mathbf{X}}_{\mathrm{Rad}}\mathbf{a}\left({\theta }_l\right),\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill {\mathbf{c}}_{\mathrm{Rad},l}& ={\left[{\mathbf{a}}^H\left({\theta }_l\right)\left\{\tilde{\mathbf{a}}\left({\varphi }_1\right){\tilde{\mathbf{a}}}^H\left({\varphi }_1\right)-{\mathbf{X}}_{\mathrm{Rad}}\right\}\mathbf{a}\left({\theta }_l\right),\cdots ,{\mathbf{a}}^H\left({\theta }_l\right)\left\{\tilde{\mathbf{a}}\left({\varphi }_U\right){\tilde{\mathbf{a}}}^H\left({\varphi }_U\right)-{\mathbf{X}}_{\mathrm{Rad}}\right\}\mathbf{a}\left({\theta }_l\right)\right]}^T.\hfill \end{array}$$

(43)

Similarly, by substituting ${\mathbf{X}}_{\mathrm{Prop}}$ into (14) for WPT operation, we have

$$\begin{array}{cc}\hfill Tr\left({\mathbf{G}}_u\left(h\right){\mathbf{X}}_{\mathrm{Prop}}\right)& =Tr\left({\mathbf{G}}_u\left(h\right){\mathbf{X}}_{\mathrm{Rad}}\right)+\sum _{k=1}^U{w}_{k}Tr\left({\mathbf{G}}_u\left(h\right)\left\{\tilde{\mathbf{a}}\left({\varphi }_k\right){\tilde{\mathbf{a}}}^H\left({\varphi }_k\right)-{\mathbf{X}}_{\mathrm{Rad}}\right\}\right)\hfill \\ \hfill & =b_{\mathrm{WPT},u}\left(h\right)+{\mathbf{c}}_{\mathrm{WPT},u}^T\left(h\right)\mathbf{w},\forall u\in \mathcal{U},\hfill \end{array}$$

(44)

where $b_{\mathrm{WPT},u}\left(h\right)\in \mathbb{R}$ and ${\mathbf{c}}_{\mathrm{WPT},u}^T\left(h\right)\in {\mathbb{R}}^{U\times 1}$ ($\forall u\in \mathcal{U}$) are functions with respect to the UAV altitude (h), expressed as

$$\begin{array}{cc}\hfill b_{\mathrm{WPT},u}\left(h\right)& =Tr\left({\mathbf{G}}_u\left(h\right){\mathbf{X}}_{\mathrm{Rad}}\right),\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill {\mathbf{c}}_{\mathrm{WPT},u}\left(h\right)& ={\left[Tr\left({\mathbf{G}}_u\left(h\right)\left\{\tilde{\mathbf{a}}\left({\varphi }_1\right){\tilde{\mathbf{a}}}^H\left({\varphi }_1\right)-{\mathbf{X}}_{\mathrm{Rad}}\right\}\right),\cdots ,Tr\left({\mathbf{G}}_u\left(h\right)\left\{\tilde{\mathbf{a}}\left({\varphi }_U\right){\tilde{\mathbf{a}}}^H\left({\varphi }_U\right)-{\mathbf{X}}_{\mathrm{Rad}}\right\}\right)\right]}^T.\hfill \end{array}$$

(46)

Accordingly, from (41) and (44), the objective function of problem ${\mathit{P}}_1$ can be rewritten as a function of $\mathbf{w}$ as follows:

$$\begin{array}{cc}\hfill & \rho L_{\mathrm{Rad}}({\mathbf{X}}_{\mathrm{Prop}},\gamma )+(1-\rho )L_{\mathrm{WPT}}({\mathbf{X}}_{\mathrm{Prop}},h)\hfill \\ \hfill & ={\displaystyle \frac{\rho }{L}}\sum _{l=1}^L{\left(\gamma \chi \left({\theta }_l\right)-b_{\mathrm{Rad},l}-{\mathbf{c}}_{\mathrm{Rad},l}^T\mathbf{w}\right)}^2+{\displaystyle \frac{1-\rho }{U}}\sum _{u=1}^U{\left(1-{\displaystyle \frac{b_{\mathrm{WPT},u}\left(h\right)+{\mathbf{c}}_{\mathrm{WPT},u}^T\left(h\right)\mathbf{w}}{P_u^{\mathrm{Des}}}}\right)}^2\hfill \\ \hfill & ={\mathbf{w}}^T\left[{\displaystyle \frac{\rho }{L}}\sum _{l=1}^L{\mathbf{c}}_{\mathrm{Rad},l}{\mathbf{c}}_{\mathrm{Rad},l}^T+{\displaystyle \frac{1-\rho }{U}}\sum _{u=1}^U{\displaystyle \frac{{\mathbf{c}}_{\mathrm{WPT},u}\left(h\right){\mathbf{c}}_{\mathrm{WPT},u}^T\left(h\right)}{{\left(P_u^{\mathrm{Des}}\right)}^2}}\right]\mathbf{w}\hfill \\ \hfill & -\left[\left(2{\displaystyle \frac{\rho }{L}}\sum _{l=1}^L\chi \left({\theta }_l\right){\mathbf{c}}_{\mathrm{Rad},l}^T\right)\gamma -2{\displaystyle \frac{\rho }{L}}\sum _{l=1}^L{b}_{\mathrm{Rad},l}{\mathbf{c}}_{\mathrm{Rad},l}^T+2{\displaystyle \frac{1-\rho }{U}}\sum _{u=1}^U\left(1-{\displaystyle \frac{b_{\mathrm{WPT},u}\left(h\right)}{P_u^{\mathrm{Des}}}}\right){\displaystyle \frac{{\mathbf{c}}_{\mathrm{WPT},u}^T\left(h\right)}{P_u^{\mathrm{Des}}}}\right]\mathbf{w}\hfill \\ \hfill & +{\displaystyle \frac{\rho }{L}}\sum _{l=1}^L{(\gamma \chi \left({\theta }_l\right)-b_{\mathrm{Rad},l})}^2+{\displaystyle \frac{1-\rho }{U}}\sum _{u=1}^U{\left(1-{\displaystyle \frac{b_{\mathrm{WPT},u}\left(h\right)}{P_u^{\mathrm{Des}}}}\right)}^2\hfill \\ \hfill & \triangleq t(\mathbf{w},\gamma ,h).\hfill \end{array}$$

(47)

As a result, the weighting factors ${\left[w_i\right]}_{i=1}^U$ for the proposed waveform, scaling factor ($\gamma$), and UAV altitude (h) can be jointly optimized by solving the following problem:

$$\begin{array}{cc}\hfill {\mathit{P}}_6:\underset{\mathbf{w},\gamma ,h}{min}& t(\mathbf{w},\gamma ,h)\hfill \end{array}$$

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& 0\le w_k\le 1,\forall k\in \mathcal{U},\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill & \sum _{k=1}^U{w}_k\le 1,\hfill \\ \hfill & H_{min}\le h\le H_{max}.\hfill \end{array}$$

(49)

Similar to the original problem, ${\mathit{P}}_1$, $t(\mathbf{w},\gamma ,h)$ in problem ${\mathit{P}}_6$ is also a non-convex, very complex, non-linear function with respect to h; therefore, it is very hard to obtain an optimal solution by directly solving it.

However, for a given $h\in [H_{min},H_{max}]$, the objective function of problem ${\mathit{P}}_6$ is convex because $t(\mathbf{w},\gamma ,h)$ is a quadratic function with respect to w; thus, we have

$$\begin{array}{c}\hfill {\displaystyle \frac{\rho }{L}}\sum _{l=1}^L{\mathbf{c}}_{\mathrm{Rad},l}{\mathbf{c}}_{\mathrm{Rad},l}^T+{\displaystyle \frac{1-\rho }{U}}\sum _{u=1}^U{\displaystyle \frac{{\mathbf{c}}_{\mathrm{WPT},u}\left(h\right){\mathbf{c}}_{\mathrm{WPT},u}^T\left(h\right)}{{\left(P_u^{\mathrm{Des}}\right)}^2}}⪰\mathbf{0},\end{array}$$

(50)

because $\rho \in [0,1]$, ${\mathbf{c}}_{\mathrm{Rad},l}{\mathbf{c}}_{\mathrm{Rad},l}^T⪰\mathbf{0}$, and ${\mathbf{c}}_{\mathrm{WPT},u}\left(h\right){\mathbf{c}}_{\mathrm{WPT},u}^T\left(h\right)⪰\mathbf{0}$. Hence, with a given h, we can obtain the optimal $\mathbf{w}$ and $\gamma$ by solving the following problem:

$$\begin{array}{cc}\hfill {\mathit{P}}_7^{\forall h}:\underset{\mathbf{w},\gamma }{min}& t(\mathbf{w},\gamma ,h)\hfill \end{array}$$

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& 0\le w_k\le 1,\forall k\in \mathcal{U},\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill & \sum _{k=1}^U{w}_k\le 1.\hfill \end{array}$$

(52)

We can numerically obtain the optimal solutions to problem ${\mathit{P}}_7^{\forall h}$ with a well-known convex optimization problem solver, such as CVX, and denote the solutions by ${\mathbf{w}}^{★}\left(h\right)$ and ${\gamma }^{★}\left(h\right)$.

Moreover, by relaxing constraints (51) and (52), the Karush–Kuhn–Tucker (KKT) optimality conditions of problem ${\mathit{P}}_7^{\forall h}$ are obtained as follows:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial t(\mathbf{w},\gamma ,h)}{\partial \mathbf{w}}}& =2\left[{\displaystyle \frac{\rho }{L}}\sum _{l=1}^L{\mathbf{c}}_{\mathrm{Rad},l}{\mathbf{c}}_{\mathrm{Rad},l}^T+{\displaystyle \frac{1-\rho }{U}}\sum _{u=1}^U{\displaystyle \frac{{\mathbf{c}}_{\mathrm{WPT},u}\left(h\right){\mathbf{c}}_{\mathrm{WPT},u}^T\left(h\right)}{{\left(P_u^{\mathrm{Des}}\right)}^2}}\right]\mathbf{w}\hfill \\ \hfill & -\left(2{\displaystyle \frac{\rho }{L}}\sum _{l=1}^L\chi \left({\theta }_l\right){\mathbf{c}}_{\mathrm{Rad},l}^T\right)\gamma +2{\displaystyle \frac{\rho }{L}}\sum _{l=1}^L{b}_{\mathrm{Rad},l}{\mathbf{c}}_{\mathrm{Rad},l}^T)\hfill \\ \hfill & -2{\displaystyle \frac{1-\rho }{U}}\sum _{u=1}^U\left(1-{\displaystyle \frac{b_{\mathrm{WPT},u}\left(h\right)}{P_u^{\mathrm{Des}}}}\right){\displaystyle \frac{{\mathbf{c}}_{\mathrm{WPT},u}^T\left(h\right)}{P_u^{\mathrm{Des}}}}\hfill \end{array}$$

(53)

$$\begin{array}{cc}\hfill & =0,\hfill \\ \hfill {\displaystyle \frac{\partial t(\mathbf{w},\gamma ,h)}{\partial \gamma }}& =-\left(2{\displaystyle \frac{\rho }{L}}\sum _{l=1}^L\chi \left({\theta }_l\right){\mathbf{c}}_{\mathrm{Rad},l}^T\right)\mathbf{w}+2{\displaystyle \frac{\rho }{L}}\sum _{l=1}^L(\gamma \chi \left({\theta }_l\right)-b_{\mathrm{Rad},l})\chi \left({\theta }_l\right)\hfill \\ & =0.\hfill \end{array}$$

(54)

From (53) and (54), we obtain the optimal ${\mathbf{w}}^{★}\left(h\right)$ and ${\gamma }^{★}\left(h\right)$ with a closed-form solution, which are expressed as

$$\begin{array}{cc}\hfill {\mathbf{w}}^{★}\left(h\right)& =\left[{\displaystyle \frac{\rho }{L}}\sum _{l=1}^L{\mathbf{c}}_{\mathrm{Rad},l}{\mathbf{c}}_{\mathrm{Rad},l}^T-{\displaystyle \frac{\rho }{L}}{\displaystyle \frac{\left({\sum }_{l=1}^L\chi \left({\theta }_l\right){\mathbf{c}}_{\mathrm{Rad},l}\right)\left({\sum }_{l=1}^L\chi \left({\theta }_l\right){\mathbf{c}}_{\mathrm{Rad},l}^T\right)}{{\sum }_{l=1}^L{\chi }^2\left({\theta }_l\right)}}\right.\hfill \\ \hfill & +{\left.{\displaystyle \frac{1-\rho }{U}}\sum _{u=1}^U{\displaystyle \frac{{\mathbf{c}}_{\mathrm{WPT},u}\left(h\right){\mathbf{c}}_{\mathrm{WPT},u}^T\left(h\right)}{{\left(P_u^{\mathrm{Des}}\right)}^2}}\right]}^{-1}\hfill \\ \hfill & \times \left[{\displaystyle \frac{\rho }{L}}{\displaystyle \frac{\left({\sum }_{l=1}^L{b}_{\mathrm{Rad},l}\chi \left({\theta }_l\right)\right)\left({\sum }_{l=1}^L\chi \left({\theta }_l\right){\mathbf{c}}_{\mathrm{Rad},l}\right)}{{\sum }_{l=1}^L{\chi }^2\left({\theta }_l\right)}}-{\displaystyle \frac{\rho }{L}}\sum _{l=1}^L{b}_{\mathrm{Rad},l}{\mathbf{c}}_{\mathrm{Rad},l}\right.\hfill \\ & +\left.{\displaystyle \frac{1-\rho }{U}}\sum _{u=1}^U\left(1-{\displaystyle \frac{b_{\mathrm{WPT},u}\left(h\right)}{P_u^{\mathrm{Des}}}}\right){\displaystyle \frac{{\mathbf{c}}_{\mathrm{WPT},u}^T\left(h\right)}{P_u^{\mathrm{Des}}}}\right],\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill {\gamma }^{★}\left(h\right)& ={\displaystyle \frac{\left({\sum }_{l=1}^L\chi \left({\theta }_l\right){\mathbf{c}}_{\mathrm{Rad},l}^T\right){\mathbf{w}}^{★}\left(h\right)+{\sum }_{l=1}^L{b}_{\mathrm{Rad},l}\chi \left({\theta }_l\right)}{{\sum }_{l=1}^L{\chi }^2\left({\theta }_l\right)}}.\hfill \end{array}$$

(56)

Here, if ${\mathbf{w}}^{★}\left(h\right)$ in (55) satisfies constraints (51) and (52) such that

$$\begin{array}{cc}\hfill & 0\le w_k^{★}\left(h\right)\le 1,\forall k\in \mathcal{U},\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill & \sum _{k=1}^U{w}_k^{★}\left(h\right)\le 1,\hfill \end{array}$$

(58)

then ${\mathbf{w}}^{★}\left(h\right)$ and ${\gamma }^{★}\left(h\right)$ are the optimal solutions to problem ${\mathit{P}}_7^{\forall h}$. Otherwise, we numerically obtain ${\mathbf{w}}^{★}\left(h\right)$ and ${\gamma }^{★}\left(h\right)$ by solving problem ${\mathit{P}}_7^{\forall h}$ when ${\mathbf{w}}^{★}\left(h\right)$ in (55) does not meet (57) and (58). Note that (55) and (56) may significantly reduce the computational complexity due to the derivation of a closed-form solution.

Similar to the procedure in Algorithm 1 in Section 3.2, by finding the optimal UAV altitude (h) based on the PSO algorithm for the one-dimensional problem, we developed a method to optimize the proposed waveform, detailed as follows: In the PSO algorithm, the m-th particle has altitude $p_m^{\left(j\right)}$ and velocity $v_m^{\left(j\right)}$ at the j-th iteration. For a given altitude $p_m^{\left(j\right)}$, the optimal weighting factor vector ${\mathbf{w}}^{★}\left(p_m^{\left(j\right)}\right)$ and scaling factor ${\gamma }^{★}\left(p_m^{\left(j\right)}\right)$ can be obtained either by applying (55) and (56) or by numerically solving problem ${\mathit{P}}_7^{\forall h}$. If the obtained ${\mathbf{w}}^{★}\left(p_m^{\left(j\right)}\right)$ from (55) satisfies constraints (57) and (58), then we use ${\mathbf{w}}^{★}\left(p_m^{\left(j\right)}\right)$ and ${\gamma }^{★}\left(p_m^{\left(j\right)}\right)$ obtained from (55) and (56) in order to further reduce the computational complexity. Otherwise, we must obtain ${\mathbf{w}}^{★}\left(p_m^{\left(j\right)}\right)$ and ${\gamma }^{★}\left(p_m^{\left(j\right)}\right)$ by solving problem ${\mathit{P}}_7^{\forall h}$. Here, the fitness function, denoted by $f_{\mathrm{Fit}}\left(p_m^{\left(j\right)}\right)$, is given by

$$\begin{array}{c}\hfill f_{\mathrm{Fit}}^{\mathrm{Prop}}\left(p_m^{\left(j\right)}\right)=t\left({\mathbf{w}}^{★}\left(p_m^{\left(j\right)}\right),{\gamma }^{★}\left(p_m^{\left(j\right)}\right),p_m^{\left(j\right)}\right).\end{array}$$

(59)

Then, the altitude $p_m^{\left(j\right)}$ and velocity $v_m^{\left(j\right)}$ are updated by minimizing $f_{\mathrm{Fit}}^{\mathrm{Prop}}\left(p_m^{\left(j\right)}\right)$. Accordingly, we can obtain the optimal altitude, $h_{\mathrm{Prop}}^{★}=g^{best}$, based on the PSO algorithm. Finally, we obtain the optimal weighting factor vector (${\mathbf{w}}^{★}\left(h_{\mathrm{Prop}}^{★}\right)$) and scaling factor (${\gamma }^{★}\left(h_{\mathrm{Prop}}^{★}\right)$) and then determine the proposed waveform, ${\mathbf{X}}_{\mathrm{Prop}}^{★}$, using (37). The detailed procedure of the proposed method is described in Algorithm 2.

Algorithm 2 Algorithm for proposed efficient design method

1: Input: A sampled angle grid ${\left[{\theta }_l\right]}_{l=1}^L$ and a desired transmit beam pattern $\chi \left({\theta }_l\right)$,             parameters on the channel statistics ($a_1,b_1,a_2,b_2,a_3,b_3$), minimum and             maximum altitude of $H_{min}$ and $H_{max}$, sum power constraint $P_T$,             all EHDs’ horizontal distance, LoS component, and desired power,             i.e., ${\left[r_u,{\mathbf{a}}_u,P_u^{\mathrm{Des}}>0\right]}_{u=1}^U$, coefficients for PSO algorithm ($\sigma ,z_1,z_2$) 2: Obtain ${\mathbf{X}}_{\mathrm{Rad}}$ by solving problem ${\mathit{P}}_5$ 3: Compute $b_{\mathrm{Rad},l}$ and ${\mathbf{c}}_{\mathrm{Rad},l}$ ($\forall l\in \mathcal{L}$) from (42) and (43), respectively 4: for Particles $m=1,2,\cdots ,M$ do 5:       Initialize $p_m^{\left(0\right)}=H_{min}+r_m^{\left(0\right)}(H_{max}-H_{min})$ and $v_m^{\left(0\right)}=0$ 6:       Initialize $p_m^{best}$←$p_m^{\left(0\right)}$ 7:       Compute $b_{\mathrm{WPT},u}\left(p_m^{best}\right)$ and ${\mathbf{c}}_{\mathrm{WPT},u}\left(p_m^{best}\right)$ ($\forall u\in \mathcal{U}$) from (45) and (46), respectively 8:       Compute ${\mathbf{w}}^{★}\left(p_m^{best}\right)$ from (55) 9:       if ${\mathbf{w}}^{★}\left(p_m^{best}\right)$ satisfies (57) and (58) then 10:           Compute ${\gamma }^{★}\left(p_m^{best}\right)$ from (56) 11:      else 12:           Solve problem ${\mathit{P}}_7^{\forall h}$ with a given $h=p_m^{best}$ and obtain ${\mathbf{w}}^{★}\left(p_m^{best}\right)$ and ${\gamma }^{★}\left(p_m^{best}\right)$ 13:      end if 14:      Compute the fitness function value $f_{\mathrm{Fit}}^{\mathrm{Prop}}\left(p_m^{best}\right)$ from (59) 15: end for 16: Initialize $g^{best}$←$arg\underset{{\left[p_m^{best}\right]}_{m=1}^M}{\min }{f}_{\mathrm{Fit}}^{\mathrm{Prop}}\left(p_m^{best}\right)$ 17: for Iterations $j=1,2,\cdots ,J$ do 18:       for Particles $m=1,2,\cdots ,M$ do 19:           Update $v_m^{\left(j\right)}$ and $p_m^{\left(j\right)}$ from (27) and (26), respectively. 20:           if $p_m^{\left(j\right)}<H_{min}$ then 21:              $p_m^{\left(j\right)}$←$H_{min}$ 22:           else if $p_m^{\left(j\right)}>H_{max}$ then 23:              $p_m^{\left(j\right)}$←$H_{max}$ 24:           end if 25:           Compute $b_{\mathrm{WPT},u}\left(p_m^{\left(j\right)}\right)$ and ${\mathbf{c}}_{\mathrm{WPT},u}\left(p_m^{\left(j\right)}\right)$ ($\forall u\in \mathcal{U}$) from (45) and (46), respectively 26:           Compute ${\mathbf{w}}^{★}\left(p_m^{\left(j\right)}\right)$ from (55) 27:           if ${\mathbf{w}}^{★}\left(p_m^{\left(j\right)}\right)$ satisfies (57) and (58) then 28:              Compute ${\gamma }^{★}\left(p_m^{\left(j\right)}\right)$ from (56) 29:           else 30:              Solve problem ${\mathit{P}}_7^{\forall h}$ with a given $h=p_m^{\left(j\right)}$ and obtain ${\mathbf{w}}^{★}\left(p_m^{\left(j\right)}\right)$ and ${\gamma }^{★}\left(p_m^{\left(j\right)}\right)$ 31:           end if 32:           Compute the fitness function value $f_{\mathrm{Fit}}^{\mathrm{Prop}}\left(p_m^{\left(j\right)}\right)$ from (59) 33:           if $f_{\mathrm{Fit}}^{\mathrm{Prop}}\left(p_m^{\left(j\right)}\right)<f_{\mathrm{Fit}}^{\mathrm{Prop}}\left(p_m^{best}\right)$ then 34:              $p_m^{best}$←$p_m^{\left(j\right)}$ 35:           end if 36:       end for 37:       Update $g^{best}$←$arg\underset{{\left[p_m^{best}\right]}_{m=1}^M}{\min }{f}_{\mathrm{Fit}}^{\mathrm{Prop}}\left(p_m^{best}\right)$ 38: end for 39: Obtain an optimal altitude $h_{\mathrm{Prop}}^{★}=g^{best}$ 40: Repeat Line 25 – 31 with $p_m^{\left(j\right)}=h_{\mathrm{Prop}}^{★}$, and  then obtain ${\mathbf{w}}^{★}\left(h_{\mathrm{Prop}}^{★}\right)$ and ${\gamma }^{★}\left(h_{\mathrm{Prop}}^{★}\right)$ 41: Obtain the proposed waveform ${\mathbf{X}}_{\mathrm{Prop}}^{★}$ from (37) 42: Output: Proposed waveform ${\mathbf{X}}_{\mathrm{Prop}}^{★}$, scaling factor ${\gamma }_{\mathrm{Prop}}^{★}$, and altitude $h_{\mathrm{Prop}}^{★}$

5. Numerical Results

In this section, we present the numerical evaluation of the proposed joint design methods. To obtain the numerical results reported in this section, we used the simulation settings listed in

Table 1. Simulation settings.

Description	Symbol	Value
Transmit sum power	$P_T$	1 W
Path loss at reference distance	${\beta }_0$	1
Path loss exponents of ground links	${\alpha }_0$	3.5
Path loss exponents of aerial links	${\alpha }_{{\displaystyle \frac{\pi }{2}}}$	2
LoS probability-related parameters	$a_1$	9.61
when considering the urban environment	$b_1$	0.16
Rician factor-related parameters	${\mathcal{K}}_0^{\mathrm{dB}}$	5 dB
	${\mathcal{K}}_{{\displaystyle \frac{\pi }{2}}}^{\mathrm{dB}}$	15 dB
Minimum altitude of UAV	$H_{min}$	50 m
Maximum altitude of UAV	$H_{max}$	250 m
Number of antennas of UAV	N	8
Coefficients in PSO algorithm	$\sigma$	0.7
	$z_1$	2
	$z_2$	2
(Number of particles)	M	8
(Number of iterations)	J	8

, similar to [5,6,7,8,9,10,11,32,33]. In addition, the carrier wavelength and the space between adjacent antennas were set to ${\displaystyle \frac{d_a}{{\lambda }_c}}=0.5$ to construct an array steering vector in (4). For the desired MIMO radar beam pattern, the sample angle grids ${\left[{\theta }_l\right]}_{l=1}^L$ were set to uniform samples in the range of $[-{90}^{\circ },{90}^{\circ }]$ with a resolution of $0.1^{\circ }$, and the desired beam pattern was formed by utilizing the dominant peak directions with its beam width [23,24]. By denoting ${\left[{\theta }_q^{\mathrm{Tar}}\right]}_{q=1}^Q$ and $\Delta$ as the dominant peak directions and beam width of the desired beam pattern ($\chi \left({\theta }_l\right)$), respectively, the latter can be formed as follows:

$$\begin{array}{c}\hfill \chi \left({\theta }_l\right)=\left\{\begin{array}{cc}1,\hfill & {\theta }_q^{\mathrm{Tar}}-{\displaystyle \frac{\Delta }{2}}\le {\theta }_l\le {\theta }_q^{\mathrm{Tar}}+{\displaystyle \frac{\Delta }{2}},q=1,2,\cdots ,Q,\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.\end{array}$$

(60)

where we set the beam width to $\Delta ={10}^{\circ }$[23,24]. Moreover, we set the desired power of the u-th EHD (i.e., $P_u^{\mathrm{Des}},\forall u\in \mathcal{U}$) to 25 % of the maximum desired power, which can be obtained by solving problem ${\mathit{P}}_4^{\prime \prime }$.

To evaluate the proposed methods, we considered the following performance metrics, which are also detailed in the figures below:

• The transmit beam pattern in terms of the angular direction with the sample angle grids ${\left[{\theta }_l\right]}_{l=1}^L$.
• The objective function value of problem ${\mathit{P}}_1$ in terms of the weighting factor $\rho$.
• The trade-off region between the beam pattern MSE, i.e., $L_{\mathrm{Rad}}(\mathbf{X},\gamma )$, and the harvested power MSE, i.e., $L_{\mathrm{WPT}}(\mathbf{X},h)$.
• The execution time for computational complexity in terms of the weighting factor $\rho$.

The transmit beam pattern refers to the spatial power in the target directions achieved by each method. The objective function value indicates the quality of the solution found by that method. Hence, lower objective function values indicate better optimization performance. The trade-off region represents the range within which improving one metric results in a compromise in the other. It demonstrates the interaction and balance between radar sensing and WPT performance, highlighting the limitations in optimizing both simultaneously. The execution time of each method serves as a practical indicator of its computational complexity, with shorter execution times signifying lower complexity.

Further, for performance comparisons, we considered the following four methods:

• The optimal waveform at UAV altitude $H_{min}$, obtained by using problem ${\mathit{P}}_2^{\forall h}$ with $h=H_{min}$.
• The optimal waveform at UAV altitude $H_{max}$, obtained by using problem ${\mathit{P}}_2^{\forall h}$ with $h=H_{max}$.
• The optimal design based on the proposed algorithm: the optimal waveform with the optimal UAV altitude obtained with Algorithm 1.
• The proposed efficient design: the proposed waveform and UAV altitude obtained with Algorithm 2.

Then, we evaluated the proposed methods with three sensing targets and two EHDs, i.e., Scenario 1 ($Q=3$ and $U=2$). By referring to the benchmark scheme for the desired beam pattern [23,24], the direction angles of the sensing targets were set to $[{\theta }_1^{\mathrm{Tar}},{\theta }_2^{\mathrm{Tar}},{\theta }_3^{\mathrm{Tar}}]=[-{40}^{\circ },0^{\circ },{40}^{\circ }]$. Regarding the EHDs’ locations, we set the horizontal distances from the UAV to $r_1=50$ m and $r_2=150$ m. In addition, to evaluate the proposed methods in various environments, we considered three different scenarios for the AoD pair of the EHDs, as follows:

• Scenario 1: $[{\varphi }_1,{\varphi }_2]=[-{20}^{\circ },{20}^{\circ }]$.
• Scenario 2: $[{\varphi }_1,{\varphi }_2]=[-{60}^{\circ },{60}^{\circ }]$.
• Scenario 3: $[{\varphi }_1,{\varphi }_2]=[-{80}^{\circ },{80}^{\circ }]$.

Figure 2 depicts the transmit beam patterns obtained with various methods using $\rho =0.6$ in Scenario 1, where $[{\varphi }_1,{\varphi }_2]=[-{20}^{\circ },{20}^{\circ }]$. For performance comparison with the optimal beam pattern for radar sensing, we also considered the radar-only design, i.e., the optimal waveform (${\mathbf{X}}_{\mathrm{Rad}}$) obtained by solving problem ${\mathit{P}}_5$.

As shown in Figure 2, with the proposed optimal and efficient designs, the peak values of the main beams in the dominant directions for radar sensing, i.e., $[-{40}^{\circ },0^{\circ },{40}^{\circ }]$, decreased compared to the radar-only design, while those of the beams in the direction of the EHDs, i.e., $[-{20}^{\circ },{20}^{\circ }]$, increased. The reason lies in the transmit power directed toward the EHDs. This result shows the effect of WPT on radar sensing. On the other hand, the peak values of the beams at $\left[{20}^{\circ }\right]$ were slightly larger than those at $[-{20}^{\circ }]$. This is because the UAV transmitted more power to the second EHD owing to the different horizontal distances of the two EHDs, i.e., $r_1\ne r_2$. Moreover, it was verified that the proposed efficient design demonstrated performance close to that of the optimal design.

Figure 3a–c show comparisons of the objective function values, the trade-off regions, and the execution times, respectively, in Scenario 1, where $[{\varphi }_1,{\varphi }_2]=[-{20}^{\circ },{20}^{\circ }]$.

Figure 3a shows that the optimal design based on the proposed method in Algorithm 1 yielded the smallest objective function value over the whole range of $\rho$. The proposed efficient design outperformed the methods with UAV altitudes $H_{min}$ and $H_{max}$, achieving near-optimal performance. In addition, Figure 3b illustrates that the proposed efficient design resulted in a considerable performance improvement compared with the methods with $H_{min}$ and $H_{max}$, achieving near-optimal performance in terms of the trade-off region. It is evident that the beam pattern MSE increased as the harvested power MSE decreased, and vice versa. Further, Figure 3c shows a comparison of the execution time of each method for evaluating the computational complexity. As expected, the proposed efficient design considerably reduced the execution time compared with the optimal design, which yielded an execution time of about 1400 s, while the proposed design had an execution time within $[20,230]$ s. It was shown that the execution time of the proposed design was comparable to those of the methods with $H_{min}$ and $H_{max}$, except when $\rho =0.9$.

Next, we conducted performance evaluations in Scenario 2, where $[{\varphi }_1,{\varphi }_2]=[-{60}^{\circ },{60}^{\circ }]$. Figure 4 shows the transmit beam patterns in this scenario.

As expected, with the optimal and proposed designs, the peak values of the main beams for radar sensing were lower than those obtained with the radar-only design. The main beams at $-{40}^{\circ }$ and ${40}^{\circ }$ were broader compared with the radar-only design. This is because the UAV transmitted signals not only for radar sensing, i.e., in the $-{40}^{\circ }$ and ${40}^{\circ }$ directions, but also in the directions of the EHDs for WPT, i.e., $-{60}^{\circ }$ and ${60}^{\circ }$. Notably, the main beam width at ${40}^{\circ }$ was much greater than that at $-{40}^{\circ }$ in order to transfer more power to the second EHD at ${60}^{\circ }$ compared with the first EHD at $-{60}^{\circ }$. On the other hand, the beam pattern obtained with the proposed efficient design was relatively different from that obtained with the optimal design, especially in the ranges of $[-{90}^{\circ },-{35}^{\circ }]$ and $[{35}^{\circ },{90}^{\circ }]$.

As such, Figure 5a–c illustrates the performance evaluations in Scenario 2, where $[{\varphi }_1,{\varphi }_2]=[-{60}^{\circ },{60}^{\circ }]$.

As shown in Figure 5a,b, the proposed efficient design outperformed the methods with $H_{min}$ and $H_{max}$ in terms of the objective function value and the trade-off region, respectively. However, it displayed performance degradation compared with the optimal design. The reason is that the main beams at $-{40}^{\circ }$ and ${40}^{\circ }$ obtained with the proposed efficient design were notably sharper than those obtained with the optimal design, as shown in Figure 4. However, the performance gap decreased as the weighting factor ($\rho$) increased. Figure 5c shows that the execution time of the optimal design was about 1400 s, whereas that of the proposed efficient design was in the range of $[58,760]$, where the time increased as $\rho$ increased. As expected, it was verified that the proposed efficient design considerably outperformed the optimal design in terms of execution time.

Finally, we conducted performance evaluations in Scenario 3, where $[{\varphi }_1,{\varphi }_2]=[-{80}^{\circ },{80}^{\circ }]$. Figure 6 illustrates the transmit beam patterns in this scenario.

As shown in Figure 6, with both the proposed optimal and efficient designs, the peak values of the main beams in the dominant directions for radar sensing, i.e., $[-{40}^{\circ },0^{\circ },{40}^{\circ }]$, were comparable to those obtained with the radar-only design. On the other hand, the beam pattern in the side-lobe regions within $[-{90}^{\circ },-{59}^{\circ }]$ and $[{57}^{\circ },{90}^{\circ }]$ was enhanced to transmit more power. The reason is that the UAV transmitted signals in the directions of the EHDs for WPT, i.e., $-{80}^{\circ }$ and ${80}^{\circ }$. Specifically, the width and power of the right side-lobe region were slightly broader and higher, respectively, than those of the left side-lobe region in order to transfer more power to the second EHD at ${80}^{\circ }$ compared with the first EHD at $-{80}^{\circ }$. Notably, it was verified that the proposed efficient design achieved similar performance to the optimal design. Further, the beam patterns of both designs closely matched those obtained with the radar-only design, except for the side-lobe regions.

Also, Figure 7a–c depict the performance evaluations for Scenario 3, where $[{\varphi }_1,{\varphi }_2]=[-{80}^{\circ },{80}^{\circ }]$.

Figure 7a,b verifies that the proposed efficient design significantly outperformed the methods with $H_{min}$ and $H_{max}$, achieving near-optimal performance in terms of both the objective function value and the trade-off region, respectively. As expected, the beam pattern MSE decreased as the harvested power MSE increased, and vice versa. Compared with those in Scenario 1, where $[{\varphi }_1,{\varphi }_2]=[-{20}^{\circ },{20}^{\circ }]$, and Scenario 2, where $[{\varphi }_1,{\varphi }_2]=[-{60}^{\circ },{60}^{\circ }]$, both the objective function value and the trade-off region in Scenario 3, where $[{\varphi }_1,{\varphi }_2]=[-{80}^{\circ },{80}^{\circ }]$, showed different performance. Additionally, Figure 7c presents a comparison of the execution time of each method. As expected, the execution time of the optimal design was about 1400 s, while that of the proposed design was within $[160,770]$. Therefore, it was shown that the proposed efficient design substantially reduced the execution time compared with the optimal design.

The numerical results verify that the proposed efficient design, which has low complexity and jointly optimizes the transmit waveform and UAV altitude, yields considerable performance improvements in terms of both the beam pattern for radar sensing and the harvested power of the EHDs, as well as the execution time due to the lower computational complexity. Additionally, the performance of UAV-enabled ISWPT systems depends on the network environment, such as the angular directions of the targets, as well as the location and desired amount of harvested power of the EHDs.

6. Conclusions

In this study, we developed a method for jointly designing the waveform and UAV altitude to accurately sense multiple targets while efficiently transferring power to multiple devices simultaneously in UAV-enabled ISWPT systems. By analyzing a challenging non-convex optimization problem, we characterized the fundamental trade-off between the desired beam pattern for radar sensing and the desired harvested power of the devices for WPT operation. The method for determining the optimal waveform and altitude was developed based on the PSO algorithm, resulting in the maximum value of the desired harvested power of the EHDs. We further developed an efficient, low-complexity method for jointly designing the transmit waveform and altitude by devising a novel transmit waveform and optimizing its parameters. According to the numerical results, it was confirmed that the proposed method offers benefits not only in terms of radar sensing and WPT performance improvements but also in reducing computational complexity. The extension of the proposed design to semi-physical simulations for practical implementation remains one of our ongoing research topics.