Joint Divergence Angle of Free Space Optics (FSO) Link and UAV Trajectory Design in FSO-Based UAV-Enabled Wireless Power Transfer Relay Systems

Abstract

Free Space Optics (FSO)-based UAV-enabled wireless power transfer (WPT) relay systems have emerged as a key technology for 6G networks, efficiently providing continuous power to Internet of Things (IoT) devices even in remote areas such as disaster recovery zones, maritime regions, and military networks, while addressing the limited battery capacity of UAVs through the FSO fronthaul link. However, the harvested power at the ground devices depends on the displacement and diameter of the FSO beam spot reaching the UAV, as well as the UAV trajectory, which affects both the FSO link and the radio-frequency (RF) link simultaneously. In this paper, we propose a joint design of the divergence angle in the FSO link and the UAV trajectory, in order to maximize the power transfer efficiency. Driven by the analysis of the optimal condition for the divergence angle, we develop a hybrid BS-PSO-based method to jointly optimize them while improving optimization performance. Numerical results demonstrate that the proposed method substantially increases power transfer efficiency and improves the optimization capability.

1. Introduction

In the 6G era, radio-frequency (RF)-based wireless power transfer (WPT) plays a vital role in enabling self-sustaining Internet of Things (IoT) networks [1,2,3]. It continuously supplies power to energy-demanding devices such as wearable electronic devices, eliminating the need for frequent battery replacements or dedicated wired power connections. Recently, as unmanned aerial vehicles (UAVs) have gained significant attention in various applications due to their flexibility in deployment, high mobility, and cost efficiency, UAV-enabled WPT has been considered a highly promising solution for efficiently delivering power to IoT devices [1,2,3]. By moving flexibly in three-dimensional (3-D) space and utilizing favorable line-of-sight (LoS) channels, UAV-enabled WPT can reduce transmission distances, avoid obstacles, and reach remote locations where conventional fixed-location power transmitters are ineffective, especially in environments such as disaster recovery, and maritime and military networks that require urgent power supply.

To reap the potential benefits in UAV-enabled WPT systems, a sophisticated UAV trajectory design is required to maximize power transfer efficiency. More specifically, the optimization of one-dimensional (1-D) and two-dimensional (2-D) trajectory design at a fixed altitude was studied to maximize the power transfer efficiency, and further extended to three-dimensional (3-D) trajectory design [4,5,6,7]. As well, a directional antenna array was studied to enhance power transfer efficiency by jointly designing the UAV trajectory and the orientation of the antenna array [8,9]. Energy beamforming, which focuses beams toward devices in the desired directions, was proposed in conjunction with the joint design of the UAV altitude to maximize power transfer efficiency [10]. Furthermore, UAV-enabled integrated sensing and WPT systems were proposed, and the joint design of the transmit waveform and UAV altitude was investigated to simultaneously improve radar sensing performance and power transfer efficiency [11]. To expand WPT coverage, a multi-UAV cooperative WPT system was studied, and an effective multi-UAV trajectory design was explored to improve power transfer performance across various environments [2,3,12,13].

On the other hand, free space optics (FSO) has been considered a key technology for 6G, as it provides high capacity by utilizing the unlicensed spectrum while strengthening security, as well as having low power consumption, reducing operational costs, and avoiding interference with the RF link [14,15,16,17,18]. The FSO system was extensively studied from a communication theory perspective, encompassing channel models, transmitter/receiver structures, modulation, channel coding, and spatial/cooperative diversity techniques [18]. A comprehensive survey of acquisition, tracking, and pointing (ATP) mechanisms for FSO systems was investigated, categorizing the mechanisms based on their working principles, use cases, and implementation technologies [14]. Also, a classification framework for FSO links according to environment, coverage, LOS availability, mobility, and distance was proposed by addressing performance variations across scenarios [15]. A comprehensive tutorial on FSO links in space was presented, covering ground-to-satellite, satellite-to-ground, and inter-satellite links [19], and was further extended to 6G non-terrestrial networks (NTNs), including airborne backhaul system architectures and their use cases [16]. Moreover, UAV-assisted FSO relay systems were extensively studied for practical demonstrations, focusing on commercially available free-flying platforms, categorized by operational altitude and payload capacity, with  an analysis of design considerations [17]. More specifically, FSO-based UAV systems that exploit the FSO link as fronthaul/backhaul link were extensively studied to enhance wireless network performance by leveraging the advantages of UAVs, which enable the establishment of strong line-of-sight (LoS) links, while handling the main challenges faced by FSO such as signal attenuation and fading due to atmospheric turbulence [20,21,22,23,24]. Not only that, the FSO-based UAV-enabled WPT system, where the UAV is supplied with power through the FSO fronthaul link to provide the RF power to the ground devices, was recently introduced to deal with the limited battery capacity issue of the UAV [3]. However, UAV trajectory design that takes into account the major challenges faced by FSO link has not been thoroughly studied despite its significant impact on WPT performance.

In the FSO-based UAV-enabled WPT system, the optical transmitter sends a very narrow beam to the UAV, and an acquisition, tracking, and pointing (ATP) system is used to precisely align the beam, guaranteeing that the UAV receives it [21,23,24]. But beam misalignment can still occur due to mechanical imperfections in the ATP system and the UAV mobility, along with unavoidable mechanical vibrations [21]. In this case, the received power at the UAV depends on the displacement and diameter of the beam spot reached at the UAV, both of which are affected by the divergence angle at the FSO transmitter. Moreover, the UAV trajectory impacts not only the received power at the UAV but also the received power at EHDs. In this context, we study the joint design of the divergence angle at the FSO transmitter and the UAV trajectory to enhance the power transfer efficiency. Our main contributions of this paper are summarized as follows:

• We formulate the problem of jointly optimizing the divergence angle and UAV trajectory to maximize the minimum harvested power among all devices to ensure fairness in FSO-based UAV-enabled WPT systems.
• To address the non-convex and highly non-linear problem, we develop a Particle Swarm Optimization (PSO)-based method to solve the problem. By leveraging the analysis on the optimal condition for the divergence angle, we further devise a hybrid BS-PSO-based method to enhance optimization performance.
• Our numerical results show that the proposed joint design substantially increases the minimum harvested power, as well as having the benefit of improving optimization capability in terms of the execution time compared to the conventional algorithm.

The rest of this paper is organized as follows. In Section 2, the related system models are introduced. In Section 3, the problem of jointly optimizing the divergence angle and UAV trajectory is formulated. In Section 4, we explore the PSO-based and the hybrid BS-PSO-based method to solve the problem. Our proposed methods are evaluated in Section 5, and a brief conclusion is drawn in Section 6.

2. System Model

Our system model is illustrated in Figure 1. We consider a free space optics (FSO)-based UAV-enabled relay system for wireless power transfer (WPT). The optical base station (OBS) transfers an optical power to the UAV via an FSO link. The UAV harvests optical signals and converts them into electrical signals, then transfers RF power to charge U energy harvesting devices (EHDs) on the ground via an RF link. Maintaining stable LoS visibility between the OBS and the UAV is crucial for providing reliable and sustainable power to the EHDs on the ground, as  any disruption in the FSO link leads to an interruption in power transfer to the UAV. In this paper, it is assumed that the LoS is maintained between the OBS and the UAV. Let us denote $(0,0,h_{\mathrm{OBS}})$ as the three-dimensional coordinates of the OBS, where $h_{\mathrm{OBS}}$ means the altitude of the OBS. Then, the 3-D locations of the UAV as well as the EHDs on the ground are represented as $(x_{\mathrm{UAV}},y_{\mathrm{UAV}},h_{\mathrm{UAV}})$ and $(x_u,y_u,0)$, where $u\in \mathcal{U}\triangleq \{1,\dots ,U\}$, respectively.

2.1. FSO Link Model Between the OBS and the UAV

By denoting $P_T$ as the transmit power of the optical beam at the FSO transmitter, the expected received optical power at the UAV, denoted as $P_R$, becomes [20,21,22]

$$\begin{array}{c}\hfill P_{\mathrm{UAV}}={\eta }_T{\eta }_R{L}_{\mathrm{Att}}\mathbb{E}\left[L_{\mathrm{Tur}}\right]\mathbb{E}\left[L_{\mathrm{Point}}\right]P_T,\end{array}$$

(1)

where ${\eta }_T$ and ${\eta }_R$ mean the optical converting efficiency between the electrical and optical power at the OBS FSO transmitter and UAV FSO receiver, respectively. According to the propagation of the optical beam, $L_{\mathrm{Att}}$, $L_{\mathrm{Tur}}$, and $L_{\mathrm{Point}}$ indicate the atmospheric attenuation loss, turbulence-induced fading, and pointing and geometric loss, respectively [20,21]. Without loss of generality, it is assumed that ${\eta }_T={\eta }_R=1$ for simplicity, where any constant values can be incorporated in our design problem. The amount of harvested energy from noise signals is ignored, as it is negligible compared to the transferred energy [1,20,21]. The details on $L_{\mathrm{Att}}$, $L_{\mathrm{Tur}}$, and $L_{\mathrm{Point}}$ are as follows.

First, in atmospheric optical beam propagation, atmospheric attenuation is mainly caused by the scattering and absorption of particles that depend on the wavelength of the optical beam. By the Beer–Lambert law, $L_{\mathrm{Atm}}$ can be modeled as [23]

$$\begin{array}{c}\hfill L_{\mathrm{Atm}}=exp(-{\sigma }_{\mathrm{Atm}}{d}_{\mathrm{FSO}}^{\mathrm{km}}),\end{array}$$

(2)

where ${\sigma }_{\mathrm{Atm}}$ denotes the atmospheric attenuation factor in dB/km units, and $d_{\mathrm{FSO}}^{\mathrm{km}}$ ($=d_{\mathrm{FSO}}/1000$) is the distance between the OBS and UAV in km units, where $d_{\mathrm{FSO}}$ (with m units) is expressed as

$$\begin{array}{c}\hfill d_{\mathrm{FSO}}=\sqrt{x_{\mathrm{UAV}}^2+y_{\mathrm{UAV}}^2+{(h_{\mathrm{OBS}}-h_{\mathrm{UAV}})}^2}.\end{array}$$

(3)

In addition, the value of ${\sigma }_{\mathrm{Atm}}$ is determined by the visibility range and wavelength of the optical beam that is given by [23]

$${\sigma }_{\mathrm{Atm}}={\displaystyle \frac{3.91}{V}}{\left({\displaystyle \frac{{\lambda }_{\mathrm{FSO}}}{550}}\right)}^{-\varrho },$$

(4)

where V and ${\lambda }_{\mathrm{FSO}}$ are the visibility range in km units and the wavelength in nm units, respectively. In (4), 550 means the reference wavelength on the visibility range with nm units. $\varrho$ is the scattering coefficient related to the particle size distribution, and its value can be determined using the Kim model as  follows [20,23]:

$$\begin{array}{c}\hfill \varrho =\left\{\begin{array}{cc}1.6,\hfill & V\ge 50,\hfill \\ 1.3,\hfill & 6<V<50,\hfill \\ 0.16V+0.34,\hfill & 1<V\le 6,\hfill \\ V-0.5,\hfill & 0.5<V\le 1,\hfill \\ 0,\hfill & V\le 0.5.\hfill \end{array}\right.\end{array}$$

(5)

Here, V depends on the weather conditions such that $V\ge 50$ with clear weather and $V\le 1$ with foggy weather [20].

Second, fading in atmospheric optical beam propagation is caused by atmospheric turbulence that is primarily due to random fluctuations in air density and temperature. In this paper, we consider the Log-normal distribution to characterize weather conditions and atmospheric turbulence effects for the FSO link [23]. Then, the probability density function (PDF) of $L_{\mathrm{Tur}}$ can be modeled as

$$\begin{array}{c}\hfill f_{L_{\mathrm{Tur}}}\left(L_{\mathrm{Tur}}\right)={\displaystyle \frac{1}{2L_{\mathrm{Tur}}\sqrt{2\pi {\xi }_x^2}}}exp\left(-{\displaystyle \frac{{\left(ln\left(L_{\mathrm{Tur}}\right)+2{\xi }_x^2\right)}^2}{8{\xi }_x^2}}\right),L_{\mathrm{Tur}}>0,\end{array}$$

(6)

where ${\xi }_x^2$ indicates the log-amplitude variance, and it is expressed as

$$\begin{array}{c}\hfill {\xi }_x^2={\displaystyle \frac{1.23C_n^2{k}^{7/6}{d}_{\mathrm{FSO}}^{11/6}}{4}},\end{array}$$

(7)

where $C_n^2$ is the refractive index structure parameter that indicates the strength of atmospheric turbulence, and $k={\displaystyle \frac{2\pi }{{\lambda }_{\mathrm{FSO}}}}$ denotes the optical wave number. From (6), the expectation of $L_{\mathrm{Tur}}$ is obtained as [21]

$$\begin{array}{cc}\hfill \mathbb{E}\left[L_{\mathrm{Tur}}\right]& ={\displaystyle \frac{1}{2}}-{\displaystyle \frac{1}{2}}\mathrm{erf}\left({\displaystyle \frac{{\xi }_x}{\sqrt{\pi }}}\right)\hfill \\ \hfill & \triangleq Q_{\mathrm{Tur}}\left(d_{\mathrm{FSO}}\right),\hfill \end{array}$$

(8)

where $\mathrm{erf}(·)$ means the error function.

At last, the details on the pointing and geometric loss $L_{\mathrm{Point}}$ are as follows. An acquisition, tracking, and pointing (ATP) system is widely adopted to precisely align a very narrow optical beam in the FSO link. However, beam misalignment may still occur due to mechanical imperfections in the ATP system and the mobility of the UAV with unavoidable mechanical vibrations [21,23,24]. By denoting $\delta$ as the the radial displacement of the optical beam at the FSO receiver, which means the distance between the center of the FSO receiver’s aperture and the center of the optical beam spot arrived at the FSO receiver, the pointing and geometric loss can be modeled as [21]

$$\begin{array}{c}\hfill L_{\mathrm{Point}}=exp\left(-{\displaystyle \frac{2{\delta }^2}{d_{\mathrm{FSO}}^2{\theta }^2}}\right)\sum _{m=0}^1{\left({\displaystyle \frac{2{\delta }^2}{d_{\mathrm{FSO}}^2{\theta }^2}}\right)}^m{\displaystyle \frac{1}{{(m!)}^2}}\gamma \left(m+1,{\displaystyle \frac{2{\mathrm{\Omega }}^2}{d_{\mathrm{FSO}}^2{\theta }^2}}\right),\end{array}$$

(9)

where $\theta$ indicates the divergence angle of the optical beam at the FSO transmitter, and $\mathrm{\Omega }$ means the radius of the FSO receiver’s aperture. In (9), $\gamma (·,·)$ represents the lower incomplete Gamma function, which is given by

$$\begin{array}{c}\hfill \gamma \left(m+1,{\displaystyle \frac{2{\mathrm{\Omega }}^2}{d_{\mathrm{FSO}}^2{\theta }^2}}\right)={\int }_0^{{\displaystyle \frac{2{\mathrm{\Omega }}^2}{d_{\mathrm{FSO}}^2{\theta }^2}}}{t}^{m}exp(-t)dt.\end{array}$$

(10)

In addition, the probability of the radial displacement $\delta$ can be modeled as the Rayleigh distribution [21,23,24], and  its PDF is represented by

$$\begin{array}{c}\hfill f_{\delta }\left(\delta \right)={\displaystyle \frac{\delta }{{\sigma }_{\mathrm{Point}}^2}}exp\left(-{\displaystyle \frac{{\delta }^2}{2{\sigma }_{\mathrm{Point}}^2}}\right),\delta \ge 0,\end{array}$$

(11)

where ${\sigma }_{\mathrm{Point}}^2$ means the variance in the radial displacement. From (9) and (11), the expectation of $L_{\mathrm{Point}}$ is obtained as [21]

$$\begin{array}{cc}\hfill \mathbb{E}\left[L_{\mathrm{Point}}\right]& ={\displaystyle \frac{1}{d_{\mathrm{FSO}}^2{\theta }^2+2{\sigma }_{\mathrm{Point}}^2}}\left[{\displaystyle \frac{d_{\mathrm{FSO}}^2{\theta }^2\left(1-g_1(d_{\mathrm{FSO}},\theta )\right)}{2}}\gamma \left(1,{\displaystyle \frac{2{\mathrm{\Omega }}^2}{d_{\mathrm{FSO}}^2{\theta }^2}}\right)\right.\hfill \\ \hfill & +\left.{\displaystyle \frac{{\sigma }_{\mathrm{Point}}^2{d}_{\mathrm{FSO}}^2{\theta }^2-g_1(d_{\mathrm{FSO}},\theta )g_2(d_{\mathrm{FSO}},\theta )}{d_{\mathrm{FSO}}^2{\theta }^2+2{\sigma }_{\mathrm{Point}}^2}}\gamma \left(2,{\displaystyle \frac{2{\mathrm{\Omega }}^2}{d_{\mathrm{FSO}}^2{\theta }^2}}\right)\right],\hfill \\ \hfill & \triangleq Q_{\mathrm{Point}}(d_{\mathrm{FSO}},\theta ),\hfill \end{array}$$

(12)

where $g_1(d_{\mathrm{FSO}},\theta )$ and $g_2(d_{\mathrm{FSO}},\theta )$ are written by

$$\begin{array}{cc}\hfill g_1(d_{\mathrm{FSO}},\theta )& =exp\left(-{\displaystyle \frac{{\left(d_{\mathrm{FSO}}\theta +\mathrm{\Omega }\right)}^2\left(d_{\mathrm{FSO}}^2{\theta }^2+2{\sigma }_{\mathrm{Point}}^2\right)}{{\sigma }_{\mathrm{Point}}^2{d}_{\mathrm{FSO}}^2{\theta }^2}}\right),\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill g_2(d_{\mathrm{FSO}},\theta )& =d_{\mathrm{FSO}}^4{\theta }^4+2\mathrm{\Omega }{d}_{\mathrm{FSO}}^3{\theta }^3\hfill \\ \hfill & +\left(3{\sigma }_{\mathrm{Point}}^2+{\mathrm{\Omega }}^2\right)d_{\mathrm{FSO}}^2{\theta }^2+4\mathrm{\Omega }{\sigma }_{\mathrm{Point}}^2{d}_{\mathrm{FSO}}\theta +2{\mathrm{\Omega }}^2{\sigma }_{\mathrm{Point}}^2.\hfill \end{array}$$

(14)

Here, we notice that $\mathbb{E}\left[L_{\mathrm{Point}}\right]$ is not a monotonic function with respect to $\theta$ [21].

2.2. RF Link Model Between the UAV and the EHDs

Next, the UAV-enabled relay system transfers RF power to U EHDs by leveraging the received optical power from the OBS, i.e., $P_{\mathrm{UAV}}$. In this case, the harvested power at the u-th EHD becomes

$$\begin{array}{cc}\hfill P_u^{\mathrm{WPT}}& ={\eta }_u{10}^{-{\alpha }_u(d_u,{\psi }_u)/10}{P}_{\mathrm{UAV}},\forall u\in \mathcal{U},\hfill \end{array}$$

(15)

where ${\eta }_u$ ($\in (0,1]$) is the energy harvesting efficiency at the u-th EHD, and ${\alpha }_u(d_u,{\psi }_u)$ is the average path loss (in dB units) for the u-th EHD’s RF link, where $d_u$ and ${\psi }_u$ denote the distance and elevation angle at the u-th EHD, respectively. Without loss of generality, we assume that ${\eta }_1=\cdots ={\eta }_U=1$ for simplicity, where any constant values can be incorporated in our design problem. The details on ${\alpha }_u(d_u,{\psi }_u)$ are as follows.

For the aerial RF link, we consider the widely adopted air-to-ground (A2G) channel model, where a probability of having a line-of-sight (LoS) link for the u-th EHD depends on the elevation angle between the UAV and it [25,26,27]. Let $l_u$ and $d_u$ be the horizontal distance and distance between the UAV and u-th EHD, respectively, then $l_u$ and $d_u$ are expressed as

$$\begin{array}{cc}\hfill l_u& =\sqrt{{(x_{\mathrm{UAV}}-x_u)}^2+{(y_{\mathrm{UAV}}-y_u)}^2},\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill d_u& =\sqrt{l_u^2+h_{\mathrm{UAV}}^2}.\hfill \end{array}$$

(17)

Then, the elevation angle (in radian units) between the UAV and u-th EHD, denoted by ${\psi }_u$, is represented by

$$\begin{array}{cc}\hfill {\psi }_u& =arctan\left({\displaystyle \frac{h_{\mathrm{UAV}}}{l_u}}\right).\hfill \end{array}$$

(18)

Here, we notice that $d_u$ and ${\psi }_u$ depend on the UAV location.

Moreover, the probability of having a LoS link for the u-th EH, denoted by $P_{\mathrm{LoS}}\left({\psi }_u\right)$, is the function of its elevation angle ${\psi }_u$, and  it is modeled by [25,26,27]

$$\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}$$

(19)

where $a_1>0$ and $b_1>0$ indicate the environmentally dependent (i.e., suburban, urban, dense urban, and highrise urban) parameters, respectively. Then, the average path loss ${\alpha }_u$ depends on both $d_u$ and ${\psi }_u$, and  it is modeled as [25,26,27]

$$\begin{array}{cc}\hfill {\alpha }_u(d_u,{\psi }_u)& =20{log}_{10}\left({\displaystyle \frac{4\pi f_0{d}_u}{c_0}}\right)+P_{\mathrm{LoS}}\left({\psi }_u\right)\times {\zeta }^{\mathrm{LoS}}+\left(1-P_{\mathrm{LoS}}\left({\psi }_u\right)\right)\times {\zeta }^{\mathrm{NLoS}},\hfill \end{array}$$

(20)

where $f_0$ and $c_0$ are the carrier frequency of the RF transmission and speed of light, respectively. In (20), ${\zeta }^{\mathrm{LoS}}$ and ${\zeta }^{\mathrm{NLoS}}$ denote the environmental dependent average additional path loss for the LoS and NLoS in dB units, respectively.

3. Problem Formulation

According to (2), (8), (12), and  (15), the average harvested power at the u-th EHD in the FSO-based UAV-enabled relay system for WPT is written by

$$\begin{array}{cc}\hfill P_u^{\mathrm{WPT}}& =P_T{10}^{-{\alpha }_u(d_u,{\psi }_u)/10}exp(-{\sigma }_{\mathrm{Atm}}{d}_{\mathrm{FSO}}/1000)Q_{\mathrm{Tur}}\left(d_{\mathrm{FSO}}\right)Q_{\mathrm{Point}}(d_{\mathrm{FSO}},\theta )\hfill \\ \hfill & \triangleq P_T{Q}_u^{\mathrm{WPT}}(x_{\mathrm{UAV}},y_{\mathrm{UAV}},h_{\mathrm{UAV}},\theta ),\hfill \end{array}$$

(21)

where $d_{\mathrm{FSO}}$, ${\left[d_u\right]}_{u=1}^U$, and ${\left[{\psi }_u\right]}_{u=1}^U$ are the functions with respect to the UAV location $(x_{\mathrm{UAV}},y_{\mathrm{UAV}},h_{\mathrm{UAV}})$. Moving the UAV closer to the OBS reduces the distance between the OBS and UAV, i.e., $d_{\mathrm{FSO}}$, but increases the distances between the UAV and the EHDs, i.e., ${\left[d_u\right]}_{u=1}^U$, and vice versa. In addition, the UAV location affects not only the distance $d_u$ of the u-th EHD but also its elevation angle ${\psi }_u$ that has an impact on the LoS probability $P_{\mathrm{LoS}}\left({\psi }_u\right)$. Furthermore, the pointing and geometric loss is affected by both the UAV location and the divergence angle $\theta$ at the OBS. Therefore, the UAV location and divergence angle need to be jointly optimized by taking account of all their interactions and their impact on the harvested power at the EHDs. To guarantee fairness in the harvested power among all EHDs, our design objective is to jointly optimize the divergence angle and the UAV trajectory to maximize the minimum harvested power among all devices. Hence, we formulate our design problem as follows:

$$\begin{array}{cc}\hfill {\mathit{P}}_1:\underset{x_{\mathrm{UAV}},y_{\mathrm{UAV}},h_{\mathrm{UAV}},\theta }{maximize}& min\left(Q_1^{\mathrm{WPT}}(x_{\mathrm{UAV}},y_{\mathrm{UAV}},h_{\mathrm{UAV}},\theta ),\dots ,Q_U^{\mathrm{WPT}}(x_{\mathrm{UAV}},y_{\mathrm{UAV}},h_{\mathrm{UAV}},\theta )\right)\hfill \end{array}$$

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& x_{min}\le x_{\mathrm{UAV}}\le x_{max},\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & y_{min}\le y_{\mathrm{UAV}}\le y_{max},\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & h_{min}\le h_{\mathrm{UAV}}\le h_{max}.\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill & {\theta }_{min}\le \theta \le {\theta }_{max},\hfill \end{array}$$

(25)

where the constraints (22)–(24) mean the coverage of the UAV flight paths. Specifically, the constraints (22) and (23) with $x_{min}$, $x_{max}$, $y_{min}$, and $y_{max}$ imply the allowable horizontal coverage, and  the constraint (24) with $h_{min}$ and $h_{max}$ means the allowable altitude for the UAV in order to receive the optical beam from the OBS by avoiding obstacles or danger zones such as buildings [28]. In addition, the constraint (25) represents the divergence angle boundary, where ${\theta }_{min}$ and ${\theta }_{max}$ represent the minimum and maximum divergence angles at the FSO transmitter [21].

The objective function of the problem ${\mathit{P}}_1$ is a very complicated non-linear function as well as a non-convex one with respect to $x_{\mathrm{UAV}},y_{\mathrm{UAV}},h_{\mathrm{UAV}},\theta$. Hence, it is extremely difficult to solve the problem in general. In the next section, we explore the Particle Swarm Optimization (PSO)-based method to find an optimal solution.

4. Joint Design of Divergence Angle of FSO Link and UAV Trajectory

We first introduce the basic idea and process of the PSO algorithm. Then, we develop the PSO-based algorithm to solve our design problem. Finally, we further devise the hybrid BS-PSO-based method to improve the optimization ability by leveraging the analysis on the optimal condition for the divergence angle.

4.1. Preliminary on Particle Swarm Optimization (PSO)

Particle Swarm Optimization (PSO) is a global optimization algorithm-based on swarm intelligence, commonly used as a meta-heuristic approach [11,28,29,30,31]. A group of particles, each with a position and velocity, explores the solution space to find the optimal solution by taking into account both its own best-discovered solution and the best solution found by the swarm as a whole. The particles iteratively update their positions and velocities while balancing exploration through individual learning and exploitation through cooperation with the whole swarm. This process is guided by a fitness function that helps refine the search path, leading to the best possible solution for optimizing the fitness function. Thanks to its efficiency with easy implementation, the PSO algorithm is widely used in various applications such as neural networks, function optimization, control, and so on [11,28,29,30,31]. Specifically, it has been extensively applied in the UAV trajectory and path planning optimizations. The details of the PSO algorithm for the L-dimensional optimization problem are as follows.

A swarm of M particles, where M denotes the total number of particles, navigates the L-dimensional variable space by updating the position of each particle at each iteration as follows:

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

(26)

where ${\mathbf{d}}_m^{\left(j\right)}=\left[d_{m,1}^{\left(j\right)},d_{m,2}^{\left(j\right)},\dots ,d_{m,L}^{\left(j\right)}\right]\in {\mathbb{R}}^L$ and ${\mathbf{v}}_m^{\left(j\right)}=\left[v_{m,1}^{\left(j\right)},v_{m,2}^{\left(j\right)},\dots ,v_{m,L}^{\left(j\right)}\right]\in {\mathbb{R}}^L$ indicate the position and velocity of m-th particle at the j-th iteration, respectively, where $m\in \{1,\dots ,M\}$ and $j\in \{0,\dots ,J-1\}$, and J means the total number of iterations [11,28,29,30,31]. In this case, ${\mathbf{v}}_m^{\left(j\right)}$ implies the rate of the next movement to update its new position. Specifically, by denoting the position of the best particle obtained in the m-th particle until the current iteration (i.e., $1,2,\dots ,j$) as ${\mathbf{d}}_m^{best}=\left[d_{m,1}^{best},d_{m,2}^{best},\dots ,d_{m,L}^{best}\right]\in {\mathbb{R}}^L$, and the position of the best particle among the whole swarm until the current iteration as ${\mathbf{g}}^{best}=\left[g_1^{best},g_2^{best},\dots ,g_L^{best}\right]\in {\mathbb{R}}^L$, the PSO algorithm updates the velocity of the m-th particle at the j-th iteration as follows:

$$\begin{array}{c}\hfill v_{m,l}^{(j+1)}=\sigma ·v_{m,l}^{\left(j\right)}+z_1·r_{m,l,1}^{(j+1)}·\left(d_{m,l}^{best}-d_{m,l}^{\left(j\right)}\right)+z_2·r_{m,l,2}^{(j+1)}·\left(g_l^{best}-d_{m,l}^{\left(j\right)}\right),l=1,2,\dots ,L,\end{array}$$

(27)

where $\sigma \in (0,1)$ represents the inertia weight, which controls the impact of momentum in the previous velocity. In addition, $z_1$ is a self-cognition coefficient for controlling the ability to learn from the particle itself, while $z_2$ is a coefficient that reflects the influence of the whole swarm on the particle. $r_{m,l,1}^{\left(j\right)}$ and $r_{m,l,2}^{\left(j\right)}$ imply the uniformly distributed random numbers in the range of $[0,1]$, respectively. As in (27), each particle’s movement is guided by its local best known position, i.e., ${\mathbf{d}}_m^{best}$, as well as its global known position, i.e., ${\mathbf{g}}^{best}$, which is updated by other, better performing particles throughout all previous iterations [11,28,29,30,31].

4.2. Joint Design Based on the PSO Method

In the proposed PSO-based optimization, a swarm of particles means the four-dimensional vector (i.e., $L=4$), which is composed of the UAV location and divergence angle. Thus, the position vector of the m-th particle at the j-th iteration is represented as

$$\begin{array}{cc}\hfill {\mathbf{d}}_m^{\left(j\right)}& ={\left[x_{\mathrm{UAV},m}^{\left(j\right)},y_{\mathrm{UAV},m}^{\left(j\right)},h_{\mathrm{UAV},m}^{\left(j\right)},{\theta }_m^{\left(j\right)}\right]}^T,\hfill \\ \hfill & \triangleq {\left[d_{m,1}^{\left(j\right)},d_{m,2}^{\left(j\right)},d_{m,3}^{\left(j\right)},d_{m,4}^{\left(j\right)}\right]}^T.\hfill \end{array}$$

(28)

In (28), each element of the position vector must satisfy the constraints (22)–(25) such that

$$\begin{array}{c}\hfill d_{min,l}\le d_{m,l}^{\left(j\right)}\le d_{max,l},l=1,2,3,4,\end{array}$$

(29)

where $d_{min,l}$ and $d_{max,l}$ are the l-th element of ${\mathbf{d}}_{min}$ and ${\mathbf{d}}_{max}$, respectively, and $d_{min,l}$ and $d_{max,l}$ are given, respectively, as 

$$\begin{array}{cc}\hfill {\mathbf{d}}_{min}& ={[x_{min},y_{min},h_{min},{\theta }_{min}]}^T,\hfill \\ \hfill {\mathbf{d}}_{max}& ={[x_{max},y_{max},h_{max},{\theta }_{max}]}^T.\hfill \end{array}$$

To find an optimal solution of the problem ${\mathit{P}}_1$, the objective function of the problem becomes the fitness function of the proposed PSO algorithm. By denoting $f_{\mathrm{Fit}}\left({\mathbf{d}}_m^{\left(j\right)}\right)$ as the fitness function, it is expressed as

$$\begin{array}{c}\hfill f_{\mathrm{Fit}}\left({\mathbf{d}}_m^{\left(j\right)}\right)=min\left(Q_1^{\mathrm{WPT}}\left(d_{m,1}^{\left(j\right)},d_{m,2}^{\left(j\right)},d_{m,3}^{\left(j\right)},d_{m,4}^{\left(j\right)}\right),\dots ,Q_U^{\mathrm{WPT}}\left(d_{m,1}^{\left(j\right)},d_{m,2}^{\left(j\right)},d_{m,3}^{\left(j\right)},d_{m,4}^{\left(j\right)}\right)\right).\end{array}$$

(30)

Then, ${\mathbf{d}}_m^{\left(j\right)}$ is iteratively updated for increasing the fitness function $f_{\mathrm{Fit}}\left({\mathbf{d}}_m^{\left(j\right)}\right)$, and the PSO algorithm is terminated when $f_{\mathrm{Fit}}\left({\mathbf{d}}_m^{\left(j\right)}\right)$ does not improve for a certain number of iterations. Finally, we obtain an best solution for maximizing $f_{\mathrm{Fit}}\left({\mathbf{d}}_m^{\left(j\right)}\right)$.

The detailed procedure is described in Algorithm 1. Lines 3–7 describe the initialization of the velocity ${\mathbf{v}}_m^{\left(0\right)}$ and the position ${\mathbf{d}}_m^{\left(0\right)}$ satisfying the constraint (29), as well as the initialization of ${\mathbf{d}}_m^{best}$, for all particles. Also, Line 8 means the initialization of ${\mathbf{g}}^{best}$ according to the fitness function (30). Lines 10–22 describe the update of velocity and position while satisfying the constraint (29), as well as ${\mathbf{d}}_m^{best}$ based on the fitness function (30). Line 23 denotes the update of ${\mathbf{g}}^{best}$ to maximize the fitness function among the whole swarm. Lines 24–29 describe the termination of the PSO algorithm when the rate of change in the fitness function, i.e.,  ${\displaystyle \frac{f_{\mathrm{Fit}}\left({\mathbf{g}}^{best}\right)-f_{\mathrm{Fit}}^{best,-1}}{f_{\mathrm{Fit}}\left({\mathbf{g}}^{best}\right)}}$, does not improved for a certain number iterations ${\mathrm{Iter}}_{\mathrm{N}.\mathrm{I}.}^{max}$ with a pre-determined tolerance $ϵ$, where $f_{\mathrm{Fit}}^{best,-1}$ and $f_{\mathrm{Fit}}\left({\mathbf{g}}^{best}\right)$ represent the previous and current fitness values, and ${\mathrm{Iter}}_{\mathrm{N}.\mathrm{I}.}^{max}$ means the pre-determined maximum number of iterations for non-improved termination.

Algorithm 1 Algorithm of the PSO-based optimization.

1: Input: UAV allowable coverage $(x_{min},x_{max})$, $(y_{min},y_{max})$, and $(h_{min},h_{max})$,    divergence angle boundary $({\theta }_{min},{\theta }_{max})$,    all EHDs horizontal distance ${[x_u,y_u]}_{u=1}^U$,    coefficients for PSO algorithm ($\sigma ,z_1,z_2$), error tolerance $ϵ$,    maximum number of iterations for non-improved termination ${\mathrm{Iter}}_{\mathrm{N}.\mathrm{I}.}^{max}$ 2: Initialize $j=0$ and ${iter}_{\mathrm{N}.\mathrm{I}.}=0$ 3: for Particles $m=1,2,\dots ,M$ do 4:     Initialize $v_{m,l}^{\left(0\right)}=0$ and $d_{m,l}^{\left(0\right)}=d_{min,l}+r_{m,l,1}^{\left(0\right)}(d_{max,l}-d_{min,l})$, where $l=1,2,3,4$ 5:     Initialize $d_{m,l}^{best}=d_{m,l}^{\left(0\right)}$, where $l=1,2,3,4$ 6:     Compute the fitness function value $f_{\mathrm{Fit}}\left({\mathbf{d}}_m^{best}\right)$ from (30) 7: end for 8: Initialize ${\mathbf{g}}^{best}=arg{\displaystyle \underset{{\left[{\mathbf{d}}_m^{best}\right]}_{m=1}^M}{max}}{f}_{\mathrm{Fit}}\left({\mathbf{d}}_m^{best}\right)$, and $f_{\mathrm{Fit}}^{best,-1}=f_{\mathrm{Fit}}\left({\mathbf{g}}^{best}\right)$ 9: while${iter}_{\mathrm{N}.\mathrm{I}.}<{\mathrm{Iter}}_{\mathrm{N}.\mathrm{I}.}^{max}$do 10:     j←$j+1$ 11:     for Particles $m=1,2,\dots ,M$ do 12:         Update ${\mathbf{v}}_m^{\left(j\right)}$ and ${\mathbf{d}}_m^{\left(j\right)}$ from (27) and (26), respectively. 13:         if $d_{m,l}^{\left(j\right)}<d_{min,l}$ (where $l=1,2,3,4$) then 14:            $d_{m,l}^{\left(j\right)}$←$d_{min,l}$ 15:         else if $d_{m,l}^{\left(j\right)}>d_{max,l}$ then 16:            $d_{m,l}^{\left(j\right)}$←$d_{max,l}$ 17:         end if 18:         Compute the fitness function value $f_{\mathrm{Fit}}\left({\mathbf{d}}_m^{\left(j\right)}\right)$ from (30) 19:         if $f_{\mathrm{Fit}}\left({\mathbf{d}}_m^{\left(j\right)}\right)>f_{\mathrm{Fit}}\left({\mathbf{d}}_m^{best}\right)$ then 20:            ${\mathbf{d}}_m^{best}$←${\mathbf{d}}_m^{\left(j\right)}$ 21:         end if 22:     end for 23:     Update ${\mathbf{g}}^{best}$←$arg{\displaystyle \underset{{\left[{\mathbf{d}}_m^{best}\right]}_{m=1}^M}{max}}{f}_{\mathrm{Fit}}\left({\mathbf{d}}_m^{best}\right)$ 24:     if ${\displaystyle \frac{f_{\mathrm{Fit}}\left({\mathbf{g}}^{best}\right)-f_{\mathrm{Fit}}^{best,-1}}{f_{\mathrm{Fit}}\left({\mathbf{g}}^{best}\right)}}<ϵ$ then 25:         ${iter}_{\mathrm{N}.\mathrm{I}.}$←${iter}_{\mathrm{N}.\mathrm{I}.}+1$ 26:     else 27:         ${iter}_{\mathrm{N}.\mathrm{I}.}$← 0 28:     end if 29:      $f_{\mathrm{Fit}}^{best,-1}$← $f_{\mathrm{Fit}}\left({\mathbf{g}}^{best}\right)$ 30: end while 31: Output: Optimal solution ${[x_{\mathrm{UAV}}^{★},y_{\mathrm{UAV}}^{★},h_{\mathrm{UAV}}^{★},{\theta }^{★}]}^T={\mathbf{g}}^{best}$

4.3. Proposed Hybrid BS-PSO Method for Joint Design

In this subsection, driven by the optimal condition of $\theta$ with a given $(x_{\mathrm{UAV}},y_{\mathrm{UAV}},h_{\mathrm{UAV}})$, we devise the hybrid BS-PSO-based optimization method to enhance the optimization ability in a conventional PSO algorithm that falls into a local solution and slow convergence [28,29,30,31] as follows.

First, from (12), the upper bound of $Q_{\mathrm{Point}}(d_{\mathrm{FSO}},\theta )$ can be obtained by

$$\begin{array}{cc}\hfill & Q_{\mathrm{Point}}(d_{\mathrm{FSO}},\theta )\hfill \\ \hfill & \le {\displaystyle \frac{d_{\mathrm{FSO}}^2{\theta }^2}{2(d_{\mathrm{FSO}}^2{\theta }^2+2{\sigma }_{\mathrm{Point}}^2)}}\left[\gamma \left(1,{\displaystyle \frac{2{\mathrm{\Omega }}^2}{d_{\mathrm{FSO}}^2{\theta }^2}}\right)+{\displaystyle \frac{2{\sigma }_{\mathrm{Point}}^2}{d_{\mathrm{FSO}}^2{\theta }^2+2{\sigma }_{\mathrm{Point}}^2}}\gamma \left(2,{\displaystyle \frac{2{\mathrm{\Omega }}^2}{d_{\mathrm{FSO}}^2{\theta }^2}}\right)\right],\hfill \\ \hfill & \triangleq {\tilde{Q}}_{\mathrm{Point}}(d_{\mathrm{FSO}},\theta ),\hfill \end{array}$$

(31)

where (31) comes from the fact that $g_1(d_{\mathrm{FSO}},\theta )>0$ and $g_2(d_{\mathrm{FSO}},\theta )>0$ in (13) and (14), respectively. Then, from (31), the first partial derivative of ${\tilde{Q}}_{\mathrm{Point}}(d_{\mathrm{FSO}},\theta )$ with respect to $\theta$ is obtained by

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial {\tilde{Q}}_{\mathrm{Point}}(d_{\mathrm{FSO}},\theta )}{\partial \theta }}={\displaystyle \frac{d_{\mathrm{FSO}}^2\theta }{d_{\mathrm{FSO}}^2{\theta }^2+2{\sigma }_{\mathrm{Point}}^2}}{\widehat{Q}}_{\mathrm{Point}}(d_{\mathrm{FSO}},\theta ),\end{array}$$

(32)

where ${\widehat{Q}}_{\mathrm{Point}}(d_{\mathrm{FSO}},\theta )$ is represented as

$$\begin{array}{cc}\hfill {\widehat{Q}}_{\mathrm{Point}}(d_{\mathrm{FSO}},\theta )=& \left(1-{\displaystyle \frac{d_{\mathrm{FSO}}^2{\theta }^2}{d_{\mathrm{FSO}}^2{\theta }^2+2{\sigma }_{\mathrm{Point}}^2}}\right)\left[\gamma \left(1,{\displaystyle \frac{2{\mathrm{\Omega }}^2}{d_{\mathrm{FSO}}^2{\theta }^2}}\right)+{\displaystyle \frac{2{\sigma }_{\mathrm{Point}}^2}{d_{\mathrm{FSO}}^2{\theta }^2+2{\sigma }_{\mathrm{Point}}^2}}\gamma \left(2,{\displaystyle \frac{2{\mathrm{\Omega }}^2}{d_{\mathrm{FSO}}^2{\theta }^2}}\right)\right]\hfill \\ \hfill & -{\displaystyle \frac{\theta }{2}}\left[\left(1+{\displaystyle \frac{4{\sigma }_{\mathrm{Point}}^2{\mathrm{\Omega }}^2}{d_{\mathrm{FSO}}^2{\theta }^2(d_{\mathrm{FSO}}^2{\theta }^2+2{\sigma }_{\mathrm{Point}}^2)}}\right){\displaystyle \frac{4{\mathrm{\Omega }}^2}{d_{\mathrm{FSO}}^2{\theta }^3}}exp\left(-{\displaystyle \frac{2{\mathrm{\Omega }}^2}{d_{\mathrm{FSO}}^2{\theta }^2}}\right)\right.\hfill \\ \hfill & \left.+{\displaystyle \frac{4{\sigma }_{\mathrm{Point}}^2{d}_{\mathrm{FSO}}^2\theta }{{(d_{\mathrm{FSO}}^2{\theta }^2+2{\sigma }_{\mathrm{Point}}^2)}^2}}\gamma \left(2,{\displaystyle \frac{2{\mathrm{\Omega }}^2}{d_{\mathrm{FSO}}^2{\theta }^2}}\right)\right].\hfill \end{array}$$

(33)

Both ${\tilde{Q}}_{\mathrm{Point}}(d_{\mathrm{FSO}},\theta )$ and its first derivative are very complicated non-linear functions with respect to $\theta$, but the optimal solution of ${\theta }^{★}$ must satisfy either

$$\begin{array}{cc}\hfill & \left(i\right){\displaystyle \frac{\partial {\tilde{Q}}_{\mathrm{Point}}(d_{\mathrm{FSO}},\theta )}{\partial \theta }}{|}_{\theta ={\theta }^{★}}=0,\mathrm{or}\hfill \\ \hfill & \left(ii\right){\theta }^{★}={\theta }_{min},\mathrm{or}\hfill \\ \hfill & \left(iii\right){\theta }^{★}={\theta }_{max}.\hfill \end{array}$$

(34)

In (34), it is readily verified that $\left(i\right)$ is equivalent to ${\widehat{Q}}_{\mathrm{Point}}(d_{\mathrm{FSO}},{\theta }^{★})=0$ from (32). We notice that ${\theta }^{★}$ that satisfies $\left(i\right)$–$\left(iii\right)$ can be obtained by the well-known Bisection (BS) line-search method [32]. In this case, the upper bound of the average harvested power at the u-th EHD becomes

$$\begin{array}{cc}\hfill P_u^{\mathrm{WPT}}& <P_T{10}^{-{\alpha }_u(d_u,{\psi }_u)/10}exp(-{\sigma }_{\mathrm{Atm}}{d}_{\mathrm{FSO}}/1000)Q_{\mathrm{Tur}}\left(d_{\mathrm{FSO}}\right){\tilde{Q}}_{\mathrm{Point}}(d_{\mathrm{FSO}},\theta )\hfill \\ \hfill & \triangleq P_T{\tilde{Q}}_u^{\mathrm{WPT}}(x_{\mathrm{UAV}},y_{\mathrm{UAV}},h_{\mathrm{UAV}},\theta ),\hfill \end{array}$$

(35)

According to the optimal condition for $\theta$, we devise the hybrid BS-PSO algorithm that explores the PSO algorithm to find an optimal $(x_{\mathrm{UAV}},y_{\mathrm{UAV}},h_{\mathrm{UAV}})$, while utilizing the Bisection method to find an optimal $\theta$. In the proposed hybrid BS-PSO algorithm, a swarm of particles becomes the three-dimensional vector (i.e., $L=3$), and then the position vector of the m-th particle at the j-th iteration is represented as

$$\begin{array}{cc}\hfill {\mathbf{d}}_m^{\left(j\right)}& ={\left[x_{\mathrm{UAV},m}^{\left(j\right)},y_{\mathrm{UAV},m}^{\left(j\right)},h_{\mathrm{UAV},m}^{\left(j\right)}\right]}^T,\hfill \\ \hfill & \triangleq {\left[d_{m,1}^{\left(j\right)},d_{m,2}^{\left(j\right)},d_{m,3}^{\left(j\right)}\right]}^T.\hfill \end{array}$$

(36)

In addition, $d_{\mathrm{FSO}}$ in (3) is rewritten as

$$\begin{array}{c}\hfill d_{\mathrm{FSO}}\left({\mathbf{d}}_m^{\left(j\right)}\right)=\sqrt{{\left(d_{m,1}^{\left(j\right)}\right)}^2+{\left(d_{m,2}^{\left(j\right)}\right)}^2+{\left(h_{\mathrm{OBS}}-d_{m,3}^{\left(j\right)}\right)}^2}.\end{array}$$

(37)

The detailed procedure of the proposed hybrid BS-PSO method is described in Algorithm 2. Line 4 means the initialization ${\mathbf{v}}_m^{\left(0\right)}$, ${\mathbf{d}}_m^{\left(0\right)}$, and  ${\mathbf{d}}_m^{best}$, same as in Algorithm 1. Line 6 describes the procedure to obtain ${\theta }_m^{★}$ with a given $d_{\mathrm{FSO}}\left({\mathbf{d}}_m^{best}\right)$ based on the Bisection method, that is described in Algorithm 3. For the hybrid BS-PSO method, according to (35), the fitness function becomes

$$\begin{array}{c}\hfill f_{\mathrm{Fit}}\left({\mathbf{d}}_m^{\left(j\right)},{\theta }_m^{★,\left(j\right)}\right)=min\left({\tilde{Q}}_1^{\mathrm{WPT}}\left(d_{m,1}^{\left(j\right)},d_{m,2}^{\left(j\right)},d_{m,3}^{\left(j\right)},{\theta }_m^{★,\left(j\right)}\right),\dots ,{\tilde{Q}}_U^{\mathrm{WPT}}\left(d_{m,1}^{\left(j\right)},d_{m,2}^{\left(j\right)},d_{m,3}^{\left(j\right)},{\theta }_m^{★,\left(j\right)}\right)\right),\end{array}$$

(38)

where ${\theta }_m^{★,\left(j\right)}$ is the obtained value from the Bisection method with a given ${\mathbf{d}}_m^{\left(j\right)}$. Line 7 computes the fitness function for the initialization. Also, Line 9 describes the initialization of ${\mathbf{g}}^{best}$ as well as ${\theta }^{best}$ to maximize the fitness function (38). Lines 13–15 describe the update of ${\mathbf{d}}_m^{\left(j\right)}$ with its velocity as well as corresponding ${\theta }_m^{★,\left(j\right)}$ that is obtained from the Bisection method with a given ${\mathbf{d}}_m^{\left(j\right)}$. In addition, Lines 16–18 mean the update of ${\mathbf{d}}_m^{best}$ as well as corresponding ${\theta }_m^{★}$ for increasing the fitness function. Line 20 implies the update of both ${\mathbf{g}}^{best}$ and ${\theta }^{best}$ to maximize the fitness function among the whole swarm. Line 21 means the termination of the proposed hybrid BS-PSO method, same as in Algorithm 1.

Algorithm 2 Algorithm of the hybrid BS-PSO-based optimization.

1: Input: Same with Line 1 in Algorithm 1 2: Initialize $j=0$ and ${iter}_{\mathrm{N}.\mathrm{I}.}=0$ 3: for Particles $m=1,2,\dots ,M$ do 4:     Initialization is the same as Line 4–5 in Algorithm 1, where $l=1,2,3$ 5:     Compute $d_{\mathrm{FSO}}\left({\mathbf{d}}_m^{best}\right)$ from (37) 6:     Obtain ${\theta }_m^{★}$ from Algorithm 3 with a given $d_{\mathrm{FSO}}\left({\mathbf{d}}_m^{best}\right)$ 7:     Compute the fitness function value $f_{\mathrm{Fit}}\left({\mathbf{d}}_m^{best},{\theta }_m^{★}\right)$ from (38) 8: end for 9: Initialize $[{\mathbf{g}}^{best},{\theta }^{best}]=arg{\displaystyle \underset{{[{\mathbf{d}}_m^{best},{\theta }_m^{★}]}_{m=1}^M}{max}}{f}_{\mathrm{Fit}}\left({\mathbf{d}}_m^{best},{\theta }_m^{★}\right)$, and $f_{\mathrm{Fit}}^{best,-1}=f_{\mathrm{Fit}}\left({\mathbf{g}}^{best},{\theta }^{best}\right)$ 10: while  ${iter}_{\mathrm{N}.\mathrm{I}.}<{\mathrm{Iter}}_{\mathrm{N}.\mathrm{I}.}^{max}$do 11:     j←$j+1$ 12:     for Particles $m=1,2,\dots ,M$ do 13:         Update ${\mathbf{d}}_m^{\left(j\right)}$ same as in Line 12–17 in Algorithm 1, where $l=1,2,3$ 14:         Obtain ${\theta }_m^{★,\left(j\right)}$ from Algorithm 3 with a given $d_{\mathrm{FSO}}\left({\mathbf{d}}_m^{\left(j\right)}\right)$ 15:         Compute the fitness function value $f_{\mathrm{Fit}}\left({\mathbf{d}}_m^{\left(j\right)},{\theta }_m^{★,\left(j\right)}\right)$ from (38) 16:         if $f_{\mathrm{Fit}}\left({\mathbf{d}}_m^{\left(j\right)},{\theta }_m^{★,\left(j\right)}\right)>f_{\mathrm{Fit}}\left({\mathbf{d}}_m^{best},{\theta }_m^{★}\right)$ then 17:            ${\mathbf{d}}_m^{best}$←${\mathbf{d}}_m^{\left(j\right)}$ and ${\theta }_m^{★}$←${\theta }_m^{★,\left(j\right)}$ 18:         end if 19:     end for 20:     Update $[{\mathbf{g}}^{best},{\theta }^{best}]$← $arg{\displaystyle \underset{{[{\mathbf{d}}_m^{best},{\theta }_m^{★}]}_{m=1}^M}{max}}{f}_{\mathrm{Fit}}\left({\mathbf{d}}_m^{best},{\theta }_m^{★}\right)$ 21:     Same as in Line 24–29 in Algorithm 1 22: end while 23: Output: Optimal solution ${[x_{\mathrm{UAV}}^{★},y_{\mathrm{UAV}}^{★},h_{\mathrm{UAV}}^{★},{\theta }^{★}]}^T=[{\mathbf{g}}^{best},{\theta }^{best}]$

Algorithm 3 The Bisection line-search method to find ${\theta }_m^{★,\left(j\right)}$.

1: Input: ${\theta }_{min}$, ${\theta }_{max}$, and  ${\mathbf{d}}_m^{\left(j\right)}$ 2: ${\theta }_a={\theta }_{min},{\theta }_b={\theta }_{max}$, and then compute ${\widehat{Q}}_{\mathrm{Point}}\left(d_{\mathrm{FSO}}\left({\mathbf{d}}_m^{\left(j\right)}\right),{\theta }_a\right)$, ${\widehat{Q}}_{\mathrm{Point}}\left(d_{\mathrm{FSO}}\left({\mathbf{d}}_m^{\left(j\right)}\right),{\theta }_b\right)$ 3: while $|{\theta }_b-{\theta }_a|>ϵ$ do 4:     ${\theta }_c$←${\displaystyle \frac{{\theta }_a+{\theta }_b}{2}}$, and then compute ${\widehat{Q}}_{\mathrm{Point}}\left(d_{\mathrm{FSO}}\left({\mathbf{d}}_m^{\left(j\right)}\right),{\theta }_c\right)$ 5:     if ${\widehat{Q}}_{\mathrm{Point}}\left(d_{\mathrm{FSO}}\left({\mathbf{d}}_m^{\left(j\right)}\right),{\theta }_a\right){\widehat{Q}}_{\mathrm{Point}}\left(d_{\mathrm{FSO}}\left({\mathbf{d}}_m^{\left(j\right)}\right),{\theta }_c\right)<0$ then 6:         ${\theta }_b$←${\theta }_c$, and  ${\widehat{Q}}_{\mathrm{Point}}\left(d_{\mathrm{FSO}}\left({\mathbf{d}}_m^{\left(j\right)}\right),{\theta }_b\right)$← ${\widehat{Q}}_{\mathrm{Point}}\left(d_{\mathrm{FSO}}\left({\mathbf{d}}_m^{\left(j\right)}\right),{\theta }_c\right)$ 7:     else 8:         ${\theta }_a$←${\theta }_c$, and  ${\widehat{Q}}_{\mathrm{Point}}\left(d_{\mathrm{FSO}}\left({\mathbf{d}}_m^{\left(j\right)}\right),{\theta }_a\right)$← ${\widehat{Q}}_{\mathrm{Point}}\left(d_{\mathrm{FSO}}\left({\mathbf{d}}_m^{\left(j\right)}\right),{\theta }_c\right)$ 9:     end if 10: end while 11: ${\theta }_m^{★,\left(j\right)}$←${\displaystyle \frac{{\theta }_a+{\theta }_b}{2}}$ 12: Output: ${\theta }_m^{★,\left(j\right)}$

5. Numerical Results

In this section, we numerically evaluate our proposed joint design. This paper explores the significance of theoretical analysis and optimization in understanding and enhancing complex systems. The proposed method and numerical results validate its capability to predict scenarios and tackle challenges, offering substantial reductions in experimentation time and cost. Extending the proposed method to experimental validation for practical implementation remains one of our ongoing research topics. By referring to the previous works [10,11,19,21,25,26,27,30,31], we consider the simulation settings listed in

Table 1. Simulation settings.

Descriptions	Symbols	Values
Transmit power at the FSO transmitter [19]	$P_T$	1 W
Wavelength of an optical beam [21]	${\lambda }_{\mathrm{FSO}}$	1550 nm
Strength of atmospheric turbulence [21]	$C_n^2$	${10}^{-13}$
Radius of the FSO receiver’s aperture [21]	$\mathrm{\Omega }$	25 mm
Minimum and maximum	${\theta }_{min}$	${10}^{-7}$ rad
divergence angle [21]	${\theta }_{max}$	0.1 mrad
Carrier frequency [25,26,27]	$f_0$	2 GHz
Environmental dependent	$a_1$	9.61
LoS probability related parameters [25,26,27]	$b_1$	0.16
Environmental dependent	${\zeta }^{\mathrm{LoS}}$	1 dB
average additional path loss [25,26,27]	${\zeta }^{\mathrm{NLoS}}$	21 dB
Minimum and maximum	$h_{min}$	50 m
altitude of the UAV [10,11]	$h_{max}$	150 m
Altitude of the OBS [10,11]	$h_{\mathrm{OBS}}$	50 m
Coefficients on PSO algorithm [30,31]	$\sigma$	0.5
	$z_1$	1.5
	$z_2$	1.5
(Number of particles)	M	32
(Iteration for non-improved termination)	${\mathrm{Iter}}_{\mathrm{N}.\mathrm{I}.}^{max}$	20

.

5.1. Performance Evaluation in a Single-EHD Scenario

First, we consider a single EHD scenario, i.e., $U=1$, when the visibility range is 50 km, i.e., $V=50$ km. Figure 2 compares the received power at the UAV by varying the divergence angle for different distances between the OBS and UAV as well as the different variance in the radial displacement. For comparison, we consider that the variance in the radial displacement is ${\sigma }_{\mathrm{Point}}^2\in \{0.01,0.04\}$, and the UAV horizontal location is $(x_{\mathrm{UAV}},0)$ m with $x_{\mathrm{UAV}}\in \{500,5000\}$, respectively. Obviously, the received power at the UAV is not a monotonic function with respect to the divergence angle, and the value is maximized at $\theta =3.8\times {10}^{-5}$ rad when $({\sigma }_{\mathrm{Point}}^2,x_{\mathrm{UAV}})=(0.01,500)$, $\theta =3.8\times {10}^{-6}$ rad when $({\sigma }_{\mathrm{Point}}^2,x_{\mathrm{UAV}})=(0.01,5000)$, $\theta =6.2\times {10}^{-5}$ rad when $({\sigma }_{\mathrm{Point}}^2,x_{\mathrm{UAV}})=(0.04,500)$, and $\theta =6.2\times {10}^{-6}$ rad when $({\sigma }_{\mathrm{Point}}^2,x_{\mathrm{UAV}})=(0.04,5000)$, respectively. Figure 2 verifies that the divergence angle to maximize the received power at the UAV depends on the variance in the radial displacement ${\sigma }_{\mathrm{Point}}^2$ as well as the UAV location.

5.2. Performance Evaluation in a Multiple-EHD Scenario

Next, we evaluate our proposed methods in a multiple-EHD scenario. We consider 32 EHDs, i.e., $U=32$, which are uniformly distributed in a rectangular area on a horizontal plane with dimensions of 200 m × 200 m as illustrated in Figure 3. For performance comparison, by denoting ${\overline{x}}_{\mathrm{EHD}}$ and ${\overline{y}}_{\mathrm{EHD}}$ as the average value of ${\left[x_u\right]}_{u=1}^U$ and ${\left[y_u\right]}_{u=1}^U$, respectively, we consider the following six methods:

• Reference Method-1: $(x_{\mathrm{UAV}},y_{\mathrm{UAV}},h_{\mathrm{UAV}})=({\overline{x}}_{\mathrm{EHD}},{\overline{y}}_{\mathrm{EHD}},h_{min})$ and $\theta ={10}^{-5}$.
• Reference Method-2: $(x_{\mathrm{UAV}},y_{\mathrm{UAV}},h_{\mathrm{UAV}})=({\overline{x}}_{\mathrm{EHD}},{\overline{y}}_{\mathrm{EHD}},h_{min})$ and $\theta ={10}^{-4}$.
• Reference Method-3: $(x_{\mathrm{UAV}},y_{\mathrm{UAV}},h_{\mathrm{UAV}})=({\overline{x}}_{\mathrm{EHD}},{\overline{y}}_{\mathrm{EHD}},h_{max})$ and $\theta ={10}^{-5}$.
• Reference Method-4: $(x_{\mathrm{UAV}},y_{\mathrm{UAV}},h_{\mathrm{UAV}})=({\overline{x}}_{\mathrm{EHD}},{\overline{y}}_{\mathrm{EHD}},h_{max})$ and $\theta ={10}^{-4}$.
• Proposed Method-1: Obtained from Algorithm 1.
• Proposed Method-2: Obtained from Algorithm 2.
• Optimal Method: Optimal UAV trajectory and divergence angle obtained by 3-D exhaustive line-search method with BS line-search method.

The reference methods mean the fixed UAV location and divergence angle, while the proposed methods obtain the optimal UAV location and divergence angle from Algorithms 1 and 2, respectively. For the optimal method, we adopt the 3-D exhaustive line-search method to find an optimal solution of $(x_{\mathrm{UAV}},y_{\mathrm{UAV}},h_{\mathrm{UAV}})$ rather than the PSO method. We set the sample grid accuracy to 0.5 m for each line search.

To examine our proposed methods in various environments, we compare the minimum harvested power among all devices, by varying the variance in the radial displacement, the horizontal distance between the OBS and the center of the rectangular area (denoted as $l_{\mathrm{Area}}$), and the visibility range, respectively.

5.2.1. Performance Evaluation by Varying the Variance in the Radial Displacement

Figure 4 compares the minimum harvested power of various methods with respect to the variance in the radial displacement (i.e., ${\sigma }_{\mathrm{Point}}^2$) when $l_{\mathrm{Area}}=5$ km and $V=50$ km. As shown, the proposed methods significantly increase the minimum harvested power compared to the reference methods and achieve the optimal performance. Moreover, the harvested power decreases as the variance increases. This is because as the variance increases, the displacement of the optical beam spot reaching at the UAV increases, thereby decreasing the received power at the UAV.

As such, Figure 5a,b depict the optical divergence angle obtained from each algorithm and the corresponding execution time, respectively, by varying the variance in the radial displacement (i.e., ${\sigma }_{\mathrm{Point}}^2$) when $l_{\mathrm{Area}}=5$ km and $V=50$ km. The optimal divergence angle increases as the variance increases. The proposed methods yield similar divergence angles compared to the optimal method. Obviously, Figure 5b demonstrates that the proposed methods considerably reduce the execution time compared to the optimal method. Moreover, Proposed Method-2 yields a smaller execution time compared to Proposed Method-1.

5.2.2. Performance Evaluation by Varying the Horizontal Distance from the Center

Figure 6 compares the minimum harvested power of various methods with respect to the horizontal distance from the center (i.e., $l_{\mathrm{Area}}$) when ${\sigma }_{\mathrm{Point}}^2=0.01$ and $V=50$ km. As expected, the proposed methods considerably outperform the reference methods and achieve the optimal performance in all ranges. It is shown that the minimum harvested power of the reference methods varies according to the horizontal distance, while the harvested power of the proposed methods is barely affected by the distance. The reason is that the visibility range is $V=50$ km, resulting in atmospheric attenuation that is comparable within the range of $[1,5]$ km, while the divergence angle and UAV trajectory are optimized accordingly.

As such, Figure 7a–c illustrate the optical divergence angle, optimal UAV trajectory, and corresponding execution time obtained from each algorithm, respectively, when ${\sigma }_{\mathrm{Point}}^2=0.01$ and $V=50$ km. The optimal divergence angle decreases as the horizontal distance increases. Meanwhile, it is shown that the proposed methods achieve the optimal divergence angle as well as similar UAV trajectory over the optimal method. As expected, the proposed methods significantly reduce the execution time compared to the optimal method, and Proposed Method-2 achieves the smallest execution time in all ranges.

Next, we evaluate the minimum harvested power of various methods for different visibility ranges, i.e., $V=1$ km.

Figure 8 compares the minimum harvested power of various methods with respect to the horizontal distance from the center (i.e., $l_{\mathrm{Area}}$) when ${\sigma }_{\mathrm{Point}}^2=0.01$ and $V=1$ km. The proposed methods outperform the reference methods, and achieve the optimal performance. In particular, the minimum harvested power non-linearly decreases as the horizontal distance increases. This is because the visibility range is $V=1$ km, which lies within the range of $[1,5]$ km. Therefore, as the horizontal distance increases, atmospheric attenuation significantly increases, further highlighting the impact of reduced visibility on the signal propagation.

As such, Figure 9a–c illustrate the optical divergence angle, optimal UAV trajectory, and corresponding execution time obtained from each algorithm, respectively, when ${\sigma }_{\mathrm{Point}}^2=0.01$ and $V=1$ km. As expected, it is shown that the proposed methods achieve the near-optimal divergence angle and UAV trajectory while significantly reducing the execution time. Meanwhile, as shown in Figure 9a, the optimal divergence angle for $V=1$ km is similar to the optimal divergence angle for $V=50$ km as depicted in Figure 7a. This is because the visibility range affects only the atmospheric attenuation in the FSO link.

6. Conclusions

This paper proposed a joint design of the divergence angle and UAV trajectory for FSO-based UAV-enabled WPT relay systems. We characterized that the harvested power at the ground devices is affected by the divergence angle at the FSO transmitter as well as the UAV trajectory that impacts both the FSO link and RF link at the same time. We formulated the design problem to maximize the minimum harvested power among all devices, and developed the PSO-based method to solve the non-convex and highly non-linear problem. By investigating the optimal condition for the divergence angle, we further proposed the BS-PSO-based method to enhance the optimization ability. The proposed method was shown to significantly increase the harvested power as well as improving the optimization capability in terms of the execution time compared to the conventional algorithm.