*Article* **Analytical Blind Beamforming for a Multi-Antenna UAV Base-Station Receiver in Millimeter-Wave Bands**

**Pingchuan Liu 1,2, Kuangang Fan 2,3,\* and Yuhang Chen 2,3**


**Abstract:** Over the last decade, unmanned aerial vehicles (UAVs) with antenna arrays have usually been employed for the enhancement of wireless communication in millimeter-wave bands. They are commonly used as aerial base stations and relay platforms in order to serve multiple users. Many beamforming methods for improving communication quality based on channel estimation have been proposed. However, these methods can be resource-intensive due to the complexity of channel estimation in practice. Thus, in this paper, we formulate an MIMO blind beamforming problem at the receivers for UAV-assisted communications in which channel estimation is omitted in order to save communication resources. We introduce one analytical method, which is called the analytical constant modulus algorithm (ACMA), in order to perform blind beamforming at the UAV base station; this relies only on data received by the antenna. The feature of the constant modulus (CM) is employed to restrict the target user signals. Algebraic operations, such as singular value decomposition (SVD), are applied to separate the user signal space from other interferences. The number of users in the region served by the UAV can be detected by exploring information in the measured data. We seek solutions that are expressible as one Kronecker product structure in the signal space; then, the beamformers that correspond to each user can be successfully estimated. The simulation results show that, by using this analytically derived blind method, the system can achieve good signal recovery accuracy, a reasonable system sum rate, and acceptable complexity.

**Keywords:** UAV base station; MIMO; millimeter-wave band; blind beamforming; signal recovery

## **1. Introduction**

Over the last few decades, wider bandwidths and more robust data transformations have always been the trends of the development of wireless communication. As one new and promising technology, millimeter-wave (mmWave) techniques have provided the possibility of applying upcoming 5G technologies, such as massive multiple-input– multiple-output (MIMO) technologies, to wireless systems [1–4].

One of the highlighted application scenarios of mmWave communication is in unmanned aerial vehicle (UAVs); such scenarios are presented with a base station, communication relay platform, or other communication enhancement equipment [1–12]. For instance, Z. Xiao et al. declared the advantages and challenges of mmWave UAV base stations (UAV-BSs) in cellular networks and explored different ways of counteracting signal blockage [2]. J. Du et al. studied UAV-BSs using multi-user massive MIMO hybrid beamforming and took an energy-saving strategy into consideration due to practical constraints, such as the size and weight of UAVs [3]. In addition to energy and power constraints, the mobility and positions of ground users are also crucial factors for UAV-assisted communication systems. In Reference [4], by exploiting the mobility and position information of a UAV system and

**Citation:** Liu, P.; Fan, K.; Chen, Y. Analytical Blind Beamforming for a Multi-Antenna UAV Base-Station Receiver in Millimeter-Wave Bands. *Sensors* **2021**, *21*, 6561. https:// doi.org/10.3390/s21196561

Academic Editor: Margot Deruyck

Received: 7 September 2021 Accepted: 28 September 2021 Published: 30 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

ground users, W. Miao et al. proposed a lightweight beamforming algorithm in order to enhance the transmission performance of 5G UAV broadcasting system. In Reference [5], by using a UAV as a relay platform, F. Jiang and A. L. Swindlehurst studied the collection of messages from co-channel users on the ground and derived an algorithm for adjusting the UAV heading to maximize the approximate ergodic sum rate of the uplink channel. In some application scenarios, single-antenna UAVs are regarded as connection points or moving antennas for the enhancement of data transmission. For example, J. Ouyang et al. used a single-antenna UAV as a data transmission relay to connect access points and base stations [6]. M. Mozaffari et al. designed a novel framework for deploying a single-antenna UAV as an antenna array in order to minimize the wireless transmission time [8]. Some other recent research using single-antenna UAV-BSs has considered other scenarios [13–15].

In the research mentioned above and in other related research, usually, linear or planar antenna arrays are employed on the UAV-BSs. It is clear that adopting MIMO systems and mmWave communication in UAV-BSs has some convincing advantages. One aspect is that mmWave systems have rich spectral resources and flexible beamforming, which can meet the great demand of the data transfer rate. On the other hand, at particular positions, multi-antenna UAV-BSs can achieve high beamforming gain towards different line-of-sight (LoS) users by using mmWave communication, thus resulting in a more satisfying system throughput and communication efficiency [16–22]. However, due to some practical system constraints, such as multipath effects, blockage, and co-channel interference from different users, a UAV-BS system will sometimes have the issues of a low signal-to-noise ratio (SNR) in data acquisition and a low communication efficiency when serving multiple users.

To overcome these shortcomings and better utilize the advantages of mmWave, various beamforming techniques have been intensively investigated in UAV communication systems [16–21,23–26]. For example, Q. Song et al. studied the joint design of beamforming and power allocation to maximize the instantaneous data rate based on an efficient sub-optimal solution based on the block-coordinate descent method [23]. Some locationbased beamforming algorithms were designed to improve the secrecy outage probability performance, maximize the network throughput, and provide flexible coverage [24–26]. In Reference [16], L. Liu et al. proposed a cooperative interference cancellation strategy motivated by the Xhaul structure to mitigate the uplink multi-beam UAV communication interference. In general communication systems, LoS transmission can efficiently improve the throughput of the UAV-BS. However, in urban areas, LoS links are prone to severe deterioration due to the complex propagation environments. In Reference [17,18], reconfigurable intelligent reflecting surfaces were used to reflect the received signals in order to enhance the data transmission quality, and passive beamforming was performed to better maximize the average achievable rate and received power. Combined with machine learning and a mean field game (MFG) control scheme, L. Li et al. proposed a joint beamforming and beam-steering method in order to build a reliable and steady connection between a UAV-BS and ground users in Reference [19]. Z. Xiao et al. utilized the artificial bee colony (ABC) algorithm to find the near-optimal beamforming vector for ground users when deploying a multi-antenna UAV-BS [20]. In Reference [21], W. Zhang et al. designed a beam-training code book for moving ground users, which can be regarded as a prior work of UAV-BS beamforming.

In this research and with the most recent beamforming techniques in UAV-BS communication systems, obtaining channel state information (CSI) is usually regarded as a crucial part of prior work. The effectiveness of channel estimation/identification may greatly affect the performance of most existing methods. In practical scenarios, especially in urban areas, specific CSI is difficult to acquire due to the complex propagation environment and multi-user interference. Another factor that greatly affects the channel estimation is the mobility of transmitters and receivers. The transmitting distance can also affect the performance of channel estimation. Accurate time-varying MIMO channel estimation can be very difficult to realize [4,7,20,27]. Moreover, the positions and mobility of ground users and UAV-BSs further increase the difficulty of channel estimation [28]. This puts extra

pressure on the timeliness of a system, which causes bad beamforming performance. For this reason, we would like to formulate a blind beamforming problem for a UAV-BS that does not depend on channel estimation or CSI.

Recently, some works explored beamforming methods that needed less specific CSI. For example, X. Li et al. derived a beamforming method for different users under the assumption of only statistical CSI, i.e., the LoS component and Rician *K*-factor [29]. L. Du et al. proposed a robust pre-coder using imperfect CSI [30]. A beamforming scheme considering limited CSI was proposed with the aid of full-diversity rotation matrices [31]. However, these beamforming methods still partially rely on channel estimation, which can be time-consuming in complicated propagation environments. To tackle the thorny problem mentioned above that is caused by complicated channel estimation, in this paper, we will employ a method called the analytical constant modulus algorithm (ACMA) in order to perform blind beamforming in UAV-BS communication systems. Differently from methods that require detailed CSI, this blind method needs only the antenna output signals (the received signals) and some of their statistical information. In other words, it saves the computational costs caused by complex channel estimation. Another advantage of the blind beamforming method is that no training sequence is needed in the data transmission process, which saves the bandwidth and UAV's battery life. The task of blind beamforming is to calculate proper weight vectors for each ground user without detailed signal and channel information [32–34].

Compared with a traditional blind beamforming method called the constant modulus algorithm (CMA), the ACMA scheme can overcome the CMA's disadvantages of a lengthy iteration process and irregular convergence to local minima [32]. The analytical scheme can successfully derive a proper beamformer for the UAV-BS system without CSI. It is, in essence, a space separation method, which needs no lengthy iterations and achieves near-optimal accuracy. Using algebraic operations, such as singular value decomposition and diagonalization, we can successfully separate the user signal space (signals originally transmitted by ground users) from interference space. By properly investigating the properties of antenna outputs, this method can achieve near-optimal calculations of the weight vectors. These well-calculated weight vectors will finally formulate blind beamformers in order to assign the power of each antenna towards the corresponding ground users. At the same time, co-channel interference among different users will be efficiently compensated, and the desired signals can be successfully separated from other signals. We will evaluate this method's performance by applying it to a UAV-BS communication system in later sections. Different propagation settings are taken into consideration. The simulation results show that, through this analytical blind method, UAV-BS beamforming in an MIMO situation can provide a reasonable system sum rate while, at the same time, guaranteeing the signal reconstruction accuracy of each signal. Its complexity is comparable to that of the CMA and other classic methods under proper settings.

The main contributions of this paper are listed below:


One potential UAV-BS application scenario can be a situation in which some specific users/target signals need higher-quality communications

The rest of this paper is organized as follows: Section 2 introduces the basic models for the UAV-BS communication system. In Section 3, we formulate a blind beamforming problem. Section 4 demonstrates the way in which we solve the formulated problem. In Section 5, we evaluate the performance of our introduced method under different settings.

*Notations:* We denote scalars, vectors, and matrices by using lowercase *a*, bold lowercase *a*, and bold uppercase letters *A*, respectively. We further denote the transpose, element-wise complex conjugation, and conjugate transpose of a matrix *A* by *A* T , *A* ∗ , and *A* <sup>H</sup>, respectively. **0** denotes an all-zero-element vector or matrix of appropriate size. Moreover, ⊗ denotes a Kronecker product, whereas pinv(*A*), SVD(*A*), and rowspan(*A*) represent a pseudo-inverse, singular value decomposition, and row span of matrix *A*. Finally, (*A*)*<sup>i</sup>* denotes the *i*-th row, *A j* the *j*-th column, and *Aij* the element in row *i* and column *j* of matrix *A*, respectively.

#### **2. Basic Model**

#### *2.1. System Model*

In this paper, we consider an MIMO uplink scenario in which a multi-antenna UAV-BS serves ground users in mmWave bands, as shown in Figure 1. In this scenario, the UAV-BS is usually equipped with a half-wavelength-spaced uniform linear array (ULA) or uniform rectangular array (URA), while each ground user is equipped with a single antenna. The ground users are distributed among distinct locations, transmitting signals of the same frequency. Due to the existence of scattering and a reflective environment, a multi-path effect is also considered in this scenario. We first introduce the well-known model for describing this signal transmission process:

$$X = HS + N\_\prime \tag{1}$$

where *<sup>X</sup>* <sup>∈</sup> <sup>C</sup>*m*×*<sup>n</sup>* is the output of the *m*-element receiver on the UAV-BS over *n* symbol snapshots, *<sup>H</sup>* <sup>∈</sup> <sup>C</sup>*m*×*<sup>d</sup>* is the channel matrix between users and the UAV-BS, and *<sup>S</sup>* <sup>∈</sup> <sup>C</sup>*d*×*<sup>n</sup>* is the source signal matrix transmitted by ground users. Each *s* T *<sup>i</sup>* = (*S*)*<sup>i</sup>* <sup>∈</sup> <sup>C</sup>*<sup>n</sup>* , *i* ∈ {1, . . . , *d*} corresponds to *<sup>d</sup>* ground users. *<sup>N</sup>* <sup>∈</sup> <sup>C</sup>*m*×*<sup>n</sup>* denotes the channel noise matrix.

To fight against the impacts of the channels on the signals, beamforming is usually applied at the transmitters or receivers (sometimes both).

**Figure 1.** Application scenario of a UAV-BS serving multiple users.

In this paper, we focus on beamforming (we can also call it beam-combining) at the UAV-BS receivers, which can be modeled as

$$
\hat{\mathbf{S}} = \mathbf{W} \mathbf{X}.\tag{2}
$$

The beamforming matrix *<sup>W</sup>* <sup>∈</sup> <sup>C</sup>*d*×*<sup>m</sup>* is applied to compensate for the channels' impacts, especially for interference cancellation, power allocation, further signal recovery processes, and so on. The rows of *W*, which are denoted as *w*<sup>T</sup> *<sup>i</sup>* = (*W*)*<sup>i</sup>* , *i* ∈ {1, . . . , *d*}, correspond to the beamforming vectors for each ground user.

#### *2.2. Channel Model*

We then establish a 3D rectangular coordinate system to better illustrate the positions of the UAV-BS and ground users. The UAV-BS is located at (*x*, *y*, *huav*), where *huav* represents the UAV's flying altitude. The position of ground user *i* is (*x<sup>i</sup>* , *y<sup>i</sup>* , 0).

The channel matrix *H* consists of channel response vectors between users and the UAV-BS. In the classic 2D situation, due to the existence of multi-path components (MPCs), each channel vector *h<sup>i</sup>* can be expressed as

$$\mathfrak{M}\_{l} = \sum\_{l=1}^{L\_{l}} \kappa\_{i,l} \mathfrak{a}(m, \theta\_{i,l}),\tag{3}$$

where *κi*,*<sup>l</sup>* represents the channel gain coefficient corresponding to the *l*-th MPC of user *i*, *a*(*m*, *θi*,*<sup>l</sup>* ) is the steering vector, *m* is the number of antenna elements, *L<sup>i</sup>* is the number of existing MPCs, and *θi*,*<sup>l</sup>* is the elevation angle of arrival (AoA) of the *l*−th MPC. Due to the spatial sparsity of the angles of arrival in the mmWave channel, different MPCs will have distinct AoAs. The steering vector can then be derived as

$$\mathbf{a}(m,\theta\_{i,l}) = \left[\mathbf{1}\_{l}\mathbf{e}^{j\pi\sin\theta\_{i,l}}, \mathbf{e}^{j\pi 2\sin\theta\_{i,l}}, \dots, \mathbf{e}^{j\pi(m-1)\sin\theta\_{i,l}}\right]^\top.\tag{4}$$

Equation (4) illustrates the 2D situation, though when applying 3D beamforming by using a uniform rectangular array (URA), the steering vector will be slightly different. For a half-spaced URA, the steering vector can be derived as

$$\mathfrak{a}(m\_{\mathbf{x}\prime}m\_{\mathbf{y}\prime}\theta\_{\mathbf{i},\mathbf{l}\prime}\phi\_{\mathbf{i},\mathbf{l}}) = \mathfrak{a}\_{\mathbf{x}}(m\_{\mathbf{x}\prime}\theta\_{\mathbf{i},\mathbf{l}\prime}\phi\_{\mathbf{i},\mathbf{l}}) \otimes \mathfrak{a}\_{\mathbf{y}}(m\_{\mathbf{y}\prime}\theta\_{\mathbf{i},\mathbf{l}\prime}\phi\_{\mathbf{i},\mathbf{l}}),\tag{5}$$

where *m<sup>x</sup>* and *m<sup>y</sup>* are the URA dimensions along the *x*−axis and *y*−axis directions, and the total array elements are *m* = *m<sup>x</sup>* ∗ *my*. *φi*,*<sup>l</sup>* is the azimuth angle of the *l*−th MPC. The geometry denoted with *a<sup>x</sup>* and *a<sup>y</sup>* can further be described as

$$\mathbf{a}\_{\mathbf{x}}(m\_{\mathbf{x}\prime}\theta\_{\mathbf{i},l\prime}\phi\_{\mathbf{i},l}) = \left[\mathbf{1}\_{\prime}\mathbf{e}^{j\pi\sin\theta\_{\mathbf{i}\prime}\cos\phi\_{\mathbf{i}\prime}}, \dots, \mathbf{e}^{j\pi(m\_{\mathbf{x}}-1)\sin\theta\_{\mathbf{i}\prime}\cos\phi\_{\mathbf{i}\prime}}\right]^{\mathsf{T}}\tag{6a}$$

$$\mathbf{a}\_{\mathcal{Y}}(m\_{\mathcal{Y}}, \theta\_{i,l}, \phi\_{i,l}) = \left[\mathbf{1}, \mathbf{e}^{j\pi\sin\theta\_{i,l}\sin\phi\_{i,l}}, \dots, \mathbf{e}^{j\pi(m\_{\mathcal{Y}}-1)\sin\theta\_{i,l}\sin\phi\_{i,l}}\right]^{\mathsf{T}}.\tag{6b}$$

Generally, the MPCs between the UAV-BS and ground users are composed of the LoS part and non-LoS (NLoS) part under the condition of no blockage. For the LoS component, the elevation angle *θi*,1 and azimuth angle *φi*,1 of user *i* can be expressed as

$$\begin{cases} \theta\_{i,1} = \arctan \frac{\sqrt{(\mathbf{x} - \mathbf{x}\_i)^2 + (y - y\_i)^2}}{h\_{uav}} \\\\ \phi\_{i,1} = \arctan \frac{y - y\_i}{\mathbf{x} - \mathbf{x}\_i} \end{cases} \tag{7}$$

which are based on the UAV-BS position and user position. For the NLoS components, these angles will be some random values for which the elevation angles *θi*,*<sup>l</sup>* generally vary from –90° to 90°, while the azimuth angles *φi*,*<sup>l</sup>* vary from –180° to 180°. In the scenario of blockage, all MPCs will be NLoS components, and the channel *H* is usually described as a classic Rayleigh fading channel.

The next crucial part in determining the channel response vector *h<sup>i</sup>* is its channel gain coefficient *κi*,*<sup>l</sup>* . For the LoS component, *κi*,1 mainly depends on the path loss, which is affected by the transmitting distance *D<sup>i</sup>* and carrier frequency *f* . The channel gain coefficient for LoS can be described as

$$\kappa\_{i,1} = \frac{1}{\left(\frac{4\pi f}{c}\right) \cdot D\_i^a} = \frac{1}{\left(4\pi\lambda\right) \cdot D\_i^a} \tag{8}$$

in which *λ* is the carrier wavelength, *D<sup>i</sup>* = p (*x* − *xi*) <sup>2</sup> + (*<sup>y</sup>* − *<sup>y</sup>i*) <sup>2</sup> + *h* 2 *uav* is the linear distance between the UAV-BS and ground user *i*, and *α* is the LoS path-loss factor. For the NLoS components, we have the following channel gain:

$$\kappa\_{i,l} = \frac{\sigma\_f}{(4\pi\lambda) \cdot D\_i^{\beta}} \,\,\,\tag{9}$$

where *σ<sup>f</sup>* is a small-scale Rayleigh fading factor, and *β* is the NLoS path-loss factor. Generally, with the existence of LoS, the channel gain will be mainly dominated by the LoS part, and other NLoS parts will sometimes be considered as interferences.

In practice, due to the mobility of the UAV and ground users, some new features, such as non-stationarity and the Doppler shift effect, should also be considered in UAV communications [27,28]. Thus, the time-varying channel response vectors for a URAequipped UAV-BS should further be expressed as

$$\mathbf{h}\_{i}(t) = \sum\_{l=1}^{L\_{i}} \kappa\_{i,l}(t) \mathbf{a}(m\_{\mathbf{x}}, m\_{y}, \theta\_{i,l}(t), \phi\_{i,l}(t)) \cdot \mathbf{e}^{j\frac{2\pi}{\lambda} \int\_{t\_{0}}^{t} f\_{i,l}(t)dt},\tag{10}$$

where *f i*,*l* is the Doppler frequency. These time-varying factors in Equation (10) are, in essence, mainly determined by the status of the UAV-BS and ground users. As mentioned in Reference [27], only a few mmWave channel sounders have the ability to measure MIMO channels in time-variant environments. MIMO channels will change due to environmental variations, even when transmitters and receivers are static, not to mention in situations, such as a moving UAV-BS serving moving ground users. Due to all these difficulties, which can greatly impact the quality of the channel estimation and the latency of the system, we focus on blind methods to compensate for signal distortion caused by the channel.

#### **3. Problem Formulation**

In general, in MIMO uplink wireless communications, the basic task of beamforming is to introduce proper weight vectors *w<sup>i</sup>* in order to achieve signal recovery against interference. Generally, non-blind beamforming methods will rely on the channel information of *H* and set the channel response matrix *H* as prior knowledge that which can be obtained through channel estimation, but this consumes resources. In this paper, the problem that we aim to solve is to derive a blind beamforming scheme at the UAV-BS receivers in order to improve the communication quality of the communication systems. In other words, we aim to calculate the proper beamforming (beam-combining) matrix *W* for optimal signal recovery at the UAV receivers without detailed information about *H*. It was shown in Reference [35] that a proper *W* usually means that the mean squared error (MSE) of the recovered signals can be minimized. Thus, we formulate the blind UAV-BS beamforming problem as a signal recovery MSE minimization problem. We denote by *W X* = *S*ˆ the recovered signals; then, the problem can be expressed as finding a proper *W* that gives

$$\min \quad \frac{1}{d \cdot n} \sum\_{i=1}^{d} \sum\_{j=1}^{n} \left( \left| \mathbf{S}\_{ij} \right|^2 - \left| \mathbf{S}\_{ij} \right|^2 \right)^2 \tag{11}$$

where *d* is the total user number, and *n* is the number of snapshots. This problem illustrates that we want to minimize the signal recovery MSE for the whole system.

In the noiseless case with perfect CSI, we can directly use *W* = pinv(*H*) as a proper beamformer, resulting in

$$\mathbf{WX} = \mathbf{\hat{S}} = \mathbf{S}.\tag{12}$$

The recovered signals of the UAV-BS are exactly the original signals sent by ground users. However, in blind situations, we have no clues about what exactly the channel *H* is. Thus, we try to solve this problem in another way, which is to explore the properties of the UAV's received signal matrix *X*.

Generally, some signals in communication systems have a constant modulus (CM) feature, such as signals using BPSK or QPSK modulations. This CM feature can provide useful information for *X*. We set the user signals with the CM feature as the target signals. For the convenience of illustration, we assume that all the target user signals have unit modulus elements via an appropriate scaler. The classic constant modulus algorithm (CMA) usually regards this signal recovery problem as a constant modulus factorization problem, that is, it makes *X* = *HS* or *S* = *W X* factorizable under ideal conditions. In essence, the recovered user signals are derived from combinations of vectors in *X*, and each recovered user signal can be expressed as ˆ*s<sup>i</sup>* <sup>T</sup> = *w*<sup>T</sup> *<sup>i</sup> X*. Consequently, each ˆ*s* T *i* can be seen as a projection onto the row span of *X*. As we now have the signal CM property, for a more distinct illustration, the problem mentioned in Equation (11) can be transformed into

$$\begin{aligned} \min & \quad \frac{1}{d \cdot n} \sum\_{i=1}^{d} \sum\_{j=1}^{n} \left( 1 - \left| \hat{\mathbf{s}}\_{ij} \right|^{2} \right)^{2} \\ \text{s.t.} & \quad \hat{\mathbf{s}}\_{i}^{\mathsf{T}} = (\hat{\mathbf{S}})\_{i} \in \text{rowspan}(\mathbf{X})\_{i} \; \forall i \in \{1, \ldots, d\}. \end{aligned} \tag{13}$$

This equation guarantees that the derived beamformer explores the CM feature in order to achieve optimal beamforming for signal recovery accuracy, and simultaneously makes sure that the recovered signals belong to specific UAV-BS receivers. We will introduce an analytical blind method called the analytical constant modulus algorithm (ACMA) in order to solve this problem in the following sections.

#### **4. Solution of the Problem**

We first assume that the ground users are transmitting independent signals with enough phase richness, which means that *S* is a full rank and the CM factorization can be unique. In practice, some trivial transformations, such as phase invariance, will sometimes cause admissible non-uniqueness, but the uniqueness of recovered signals *S*ˆ can be guaranteed for *n* → ∞. Once the uniqueness is guaranteed, we can make sure that the recovered user signals at the UAV-BS site are the target signals originally sent by ground users [32]. To achieve a minimal signal recovery MSE for the whole system, we can pay attention to each user signal. Ideally, each recovered user signal should have unit modulus elements, as we assumed before, that satisfy

$$\mathbf{s}\_{i}^{\mathsf{T}} = [1, \ldots, 1]\_{\prime} \; \forall i \in \{1, \ldots, d\} \tag{14a}$$

$$\text{s.t.} \quad \hat{\mathbf{s}}\_{i}^{\mathsf{T}} = (\hat{\mathbf{S}})\_{i} \in \text{rowspan}(\mathbf{X}), \forall i \in \{1, \ldots, d\}. \tag{14b}$$

As mentioned in the previous sections, the UAV-BS's received signal matrix *X* is the only data that we can acquire in the system. For an antenna-given *<sup>X</sup>* <sup>∈</sup> <sup>C</sup>*m*×*<sup>n</sup>* , the first step in the ACMA scheme is to separate the signal space from the interference space. To divide spaces in *X*, we need to apply singular value decomposition (SVD) to *X*.

Let SVD(*X*) = *U***Σ***V* ∗ , where *<sup>U</sup>* <sup>∈</sup> <sup>C</sup>*m*×*m*, **<sup>Σ</sup>** <sup>∈</sup> <sup>R</sup>*m*×*<sup>n</sup>* , and *<sup>V</sup>* <sup>∈</sup> <sup>C</sup>*n*×*<sup>n</sup>* . By estimating the rank of *X*, we can first obtain the number of received signals. One direct way is to check how many singular values are non-zero in **Σ**; the number of non-zero values is equal to the number of transmitted signals. Remember that ground users are transmitting *d*

independent signals, and the received signal matrix *X* is, in essence, the combination of transmitted signals *s<sup>i</sup>* through channels, so there are actually *d* independent row vectors in the row space of *X*. The SVD operation gives us the singular non-zero values and their corresponding vectors, which can actually form the basis of rowspan (*X*). This SVD of *X* first gives us the number of transmitted signals; then, we extract the signal space from *X*, which is helpful for the following operations. Notice that the antenna array will sometimes receive all of the signals that appear in the system, including non-CM (non-user) ones, so the rank of *X* is not always precisely *d* in the presence of non-CM signals. We use *δ* as the total signal number; then, we have *δ* non-zero values in **Σ**. In fact, the SVD of *X* can only extract the signal space, not the target CM signal space. To obtain the CM signal space, more operations are needed.

Let *<sup>V</sup>*<sup>ˆ</sup> <sup>∈</sup> <sup>C</sup>*δ*×*<sup>n</sup>* be the submatrix of *<sup>V</sup>* that contains the singular vectors corresponding to *δ* singular non-zero values; then, we have that the rows of *V*ˆ are actually the orthogonal basis of rowspan(*X*), which can be seen as the signal space. Thus, we have

$$\mathbf{s}\_{\mathbf{i}}^{\mathsf{T}} \in \text{rowspan}(\mathbf{X}) \Leftrightarrow \mathbf{s}\_{\mathbf{i}}^{\mathsf{T}} = \omega\_{\mathbf{i}}^{\mathsf{T}} \hat{\mathbf{V}}, \quad \forall i \in \{1, \ldots, \delta\}, \tag{15}$$

where the weight vectors *w<sup>i</sup>* are acting on the signal space. Notice that, here, *s* T *i* , ∀*i* ∈ {1, . . . , *δ*} are not yet the recovered target user signals, and some non-CM signals are still in the subspace. In other words, the SVD operation of *X* extracts the *d* user signals together with *δ* − *d* interference signals. Clearly, this cannot yet meet the requirements for user signal recovery. In the following, we will separate user signals from interference signals by using the CM property.

We know that *V*ˆ = - *v*1, . . . , *v<sup>n</sup>* , where *v<sup>j</sup>* , <sup>∀</sup>*<sup>j</sup>* ∈ {1, . . . , *<sup>n</sup>*} is the *<sup>j</sup>*−th column of *<sup>V</sup>*<sup>ˆ</sup> ; then, we can rewrite Equation (14a) as

$$\begin{split} \mathfrak{s}\_{i}^{\sf T} &= [1, \ldots, 1], \forall i \in \{1, \ldots, d\} \\ &\Leftrightarrow \left[ \left| (\mathfrak{s}\_{i}^{\sf T})\_{1} \right|^{2} \ldots \left| (\mathfrak{s}\_{i}^{\sf T})\_{n} \right|^{2} \right] = [1, \ldots, 1], \forall i \in \{1, \ldots, d\} \\ &\Leftrightarrow \mathfrak{w}\_{i}^{\sf T} \mathfrak{v}\_{j} \mathfrak{v}\_{j}^{\sf H} \mathfrak{w}\_{i}^{\*} = 1, \quad \forall i \in \{1, \ldots, d\}, \forall j \in \{1, \ldots, n\}. \end{split} \tag{16}$$

We first define *P<sup>j</sup>* = *vjv* H *<sup>j</sup>* <sup>∈</sup> <sup>C</sup>*δ*×*<sup>δ</sup>* , ∀*j* ∈ {1, . . . , *n*}; then, the CM property in Equation (16) can be expressed as

$$\hat{\mathbf{s}}\_{i} = \begin{bmatrix} \boldsymbol{w}\_{i}^{\mathsf{T}} \boldsymbol{P}\_{1} \boldsymbol{w}\_{i}^{\*} \\ \vdots \\ \boldsymbol{w}\_{i}^{\mathsf{T}} \boldsymbol{P}\_{j} \boldsymbol{w}\_{i}^{\*} \\ \vdots \\ \boldsymbol{w}\_{i}^{\mathsf{T}} \boldsymbol{P}\_{h} \boldsymbol{w}\_{i}^{\*} \end{bmatrix} = \begin{bmatrix} 1 \\ \vdots \\ 1 \end{bmatrix}. \tag{17}$$

Now, in the extracted signal space, our goal is to find linearly independent beamformers *w<sup>i</sup>* for each user that can enforce this signal CM property. Once we have found all of the *w<sup>i</sup>* , we will have separated the user signal space (the target CM signal space) from the non-user signal space, and the blind beamforming problem for the UAV-BS can be solved.

For a better illustration, we can transform Equation (17) into a matrix computation form. In the following, we will temporarily drop the index *i* of *w* for clarity, i.e., *w* stands for all solutions *w<sup>i</sup>* , ∀*i* ∈ {1, . . . , *d*}. We define two operators, *vec*(·) and *mat*(·), where *vec*(·) will stack all elements of a matrix into a vector by column, and *mat*(·) is the inverse operator of *vec*(·).

Let *<sup>p</sup><sup>j</sup>* <sup>=</sup> *vec*(*Pj*) <sup>∈</sup> <sup>C</sup>*<sup>δ</sup>* 2 ; with the structure of *y* = *vec w*∗*w*<sup>T</sup> <sup>=</sup> *<sup>w</sup>* <sup>⊗</sup> *<sup>w</sup>*<sup>∗</sup> <sup>∈</sup> <sup>C</sup>*<sup>δ</sup>* 2 , we have the original CM elements in which *w*T*Pjw*<sup>∗</sup> = 1 is equal to *p* T *j y* = 1; then, Equation (17) can be expressed as

$$\mathbf{P}\mathbf{y} = \begin{bmatrix} 1 \\ \vdots \\ 1 \end{bmatrix}, \\ where \ \mathbf{P} = \begin{bmatrix} p\_1^\top \\ \vdots \\ p\_j^\top \\ \vdots \\ p\_n^\top \end{bmatrix}, \mathbf{y} = \mathbf{w} \otimes \mathbf{w}^\*. \tag{18}$$

For greater convenience of computation, we can try to use some linear algebraic methods to transform Equation (18) into a more manageable form. One reasonable choice can be Householder transformation (HT) [32,36], though other options that can achieve the same goal are also acceptable. Let

$$\mathbf{Q} = \mathbf{I} - 2 \frac{q \mathbf{q}^\*}{q^\* \mathbf{q}'} \mathbf{q} = \begin{bmatrix} 1 \\ \vdots \\ 1 \end{bmatrix} - \begin{bmatrix} \sqrt{n} \\ 0 \\ \vdots \\ 0 \end{bmatrix} . \tag{19}$$

Then, by applying this *Q* to *P*, we have separated the system into two parts:

$$\mathbf{Q}\mathbf{P}\mathbf{y} = \begin{bmatrix} \hat{\mathbf{p}}\_1 \\ \hat{\mathbf{P}} \end{bmatrix} \mathbf{y} \Leftrightarrow \begin{cases} (1) \ \hat{\mathbf{p}}\_1 \mathbf{y} = \sqrt{n} \quad \hat{\mathbf{p}}\_1 : 1 \times \boldsymbol{\delta}^2 \\ (2) \ \hat{\mathbf{P}}\mathbf{y} = \mathbf{0} \qquad \text{'} \ \hat{\mathbf{P}} : (n-1) \times \boldsymbol{\delta}^2 \end{cases} \tag{20}$$

After this transformation, we can directly focus on the second part of this system, which is easier to solve than the system *Py* = 1 . The first part corresponds to the scaling of the vectors, which will be addressed later. Now, we try to find a solution *y* that satisfies

$$
\hat{P}\mathcal{Y} = \begin{bmatrix} 0\\ \vdots\\ 0 \end{bmatrix} \prime \mathcal{Y} = \mathcal{w} \otimes \mathcal{w}^\*. \tag{21}
$$

Due to the existence of multiple ground users, there will be more than one solution to the system. We require the number of snapshots *n* to be larger than *δ* 2 ; otherwise, there will be no solutions to the system. One simple and direct way to calculate these solutions is to employ SVD for *P*ˆ because the solutions that make *Py*ˆ = **0** can be seen as the orthogonal basis of the null-space (*P*ˆ ). Thus, the singular vectors corresponding to singular zero values can be seen as the solutions. In the case of noisy channels, we retrieve the singular vectors corresponding to the smallest singular values of *P*ˆ . This operation can correctly estimate the number of user signals in the signal space.

However, as we mentioned before, each CM solution *y* corresponding to each ground user in the UAV-BS system should have the specific structure *w* ⊗ *w*<sup>∗</sup> . In common wireless UAV-BS communications, the sample size is generally much larger than the number of users served, which leads to an overdetermined system. We require *n* > *δ* 2 ; hence, the set of independent solutions *y* is not unique. Thus, we need the condition *y* = *w* ⊗ *w*<sup>∗</sup> to restrict the solution space. Simple SVD operations cannot guarantee this structure; it is not appropriate to compute the solutions *y* and hope that they have the specified structure [32]. Operations are needed to make sure that each solution can be expressed as *w* ⊗ *w*<sup>∗</sup> . We will simply introduce the extended QZ iteration method in the ACMA to solve this structural problem.

In general, the CM solution space that satisfies Equation (21) can be written as

$$\mathbf{y} = c\_1 \mathfrak{G}\_1 + c\_2 \mathfrak{G}\_2 + \dots + c\_d \mathfrak{G}\_d = \mathfrak{w} \otimes \mathfrak{w}^\*,\tag{22}$$

where {*c*1, . . . , *c* <sup>ˆ</sup>*<sup>d</sup>* } is a set of constants, {*β*1, . . . , *β* <sup>ˆ</sup>*<sup>d</sup>* } is the basis of the kernel of *<sup>P</sup>*<sup>ˆ</sup> , and <sup>ˆ</sup>*<sup>d</sup>* is the dimension of *Ker*(*P*ˆ). The linearly independent *y* corresponds to linearly independent *w*, which now leads to a linearly independent set of constants. Let SVD(*P*ˆ) = *Up***Σ***pV<sup>p</sup>* ∗ , as mentioned above, be the basis of *Ker*(*P*ˆ), which can be estimated as the singular vectors corresponding to singular zero values. The dimension of *Ker*(*P*ˆ) is equal to the number of singular zero values, which is also the number of target user signals, so we have ˆ*d* = *d* (proved in Reference [32]). Then, the last *d* rows of *V<sup>p</sup>* ∗ corresponding to singular zero values can form a basis {*β*1, . . . , *βd*}. After finding the basis of this CM solution space, once we obtain {*c*1, . . . , *cd*}, we will have separated the CM signal space from the non-CM ones. Notice that, here, we say "singular zero values" because this is under the ideal conditions, which means that there is no noise in the system. When additive noise is applied to the system, there will be a threshold for deciding when small singular values can be seen as "singular zero values". The discussion about the threshold will be introduced in the next section.

Let *B<sup>k</sup>* = *mat*(*β<sup>k</sup>* ), *k* ∈ {1, . . . , *d*}, and *Y* = *mat*(*y*); then, Equation (22) can be expressed as

$$\mathbf{Y} = \mathbf{c}\_1 \mathbf{B}\_1 + \dots + \mathbf{c}\_d \mathbf{B}\_d = \mathbf{w}^\* \cdot \mathbf{w}^\top. \tag{23}$$

A proper set of constants can make the result of Equation (23) rank 1 and positive semidefinite; hence, it can be factorized as *w*∗ · *<sup>w</sup>*T. Assume that we already have the weight vectors *w*<sup>1</sup> . . . *w<sup>d</sup>* ; then, *β<sup>i</sup>* = *w<sup>i</sup>* ⊗ *w<sup>i</sup>* ∗ , *<sup>i</sup>* ∈ {1, . . . , *<sup>d</sup>*} can form the basis of *Ker*(*P*ˆ), and each *Y<sup>i</sup>* can be factorized as

$$\mathbf{Y}\_{\mathrm{i}} = \lambda\_{\mathrm{i}1} \mathbf{w}\_{1}{}^{\*} \mathbf{w}\_{1}^{\mathsf{T}} + \lambda\_{\mathrm{i}2} \mathbf{w}\_{2}{}^{\*} \mathbf{w}\_{2}^{\mathsf{T}} + \dots + \lambda\_{\mathrm{i}d} \mathbf{w}\_{d}{}^{\*} \mathbf{w}\_{d}^{\mathsf{T}} = \mathbf{W}^{\*} \Lambda\_{\mathrm{i}} \mathbf{W}^{\mathsf{T}}.\tag{24}$$

If we could find matrix *Q* and *Z* to make *QYiZ* = **Λ***<sup>i</sup>* , then the set of constants {*c*1, . . . , *cd*} corresponding to each user signal (making the solution rank 1) could be derived from the eigenvalue matrices **Λ***<sup>i</sup>* . Consequently, the UAV-BS beamforming problem can be seen as an eigenvalue problem. To reduce the calculation complexity and save the calculation cost in the UAV-BS, it is not necessary to fully diagonalize *Y<sup>i</sup>* , as eigenvalues can also be obtained if *QYiZ* are upper-triangular matrices.

The extended QZ iteration sets the initial *Q*(0) = *Z* (0) <sup>=</sup> *<sup>I</sup>*. For the *<sup>k</sup>*−th iteration, this method aims to make

$$\begin{array}{ll}\min & ||\mathbf{Q}^{(k)}(\mathbf{Y}\_{1}\mathbf{Z}^{(k-1)})||\_{LF}^{2} + \dots + ||\mathbf{Q}^{(k)}(\mathbf{Y}\_{d}\mathbf{Z}^{(k-1)})||\_{LF}^{2} \\ \text{min} & ||(\mathbf{Q}^{(k)}\mathbf{Y}\_{1})\mathbf{Z}^{(k)}||\_{LF}^{2} + \dots + ||(\mathbf{Q}^{(k)}\mathbf{Y}\_{d})\mathbf{Z}^{(k)}||\_{LF}^{2} \end{array} \tag{25}$$

which necessitates making the lower-triangular part the minimal norm. Methods, such as HT and SVD, can be used to construct upper-triangular *R<sup>i</sup>* = *QYiZ* [32]. With the derived *Ri* , we extract the diagonal entries of each *R<sup>i</sup>* and form a coefficient matrix:

$$\mathbf{C} = \begin{bmatrix} \{\mathbf{R}\_1\}\_{11} & \dots & \{\mathbf{R}\_1\}\_{dd} \\ \vdots & & \vdots \\ \{\mathbf{R}\_d\}\_{11} & \dots & \{\mathbf{R}\_d\}\_{dd} \end{bmatrix}^{-1} \begin{array}{c} \\ \vdots \\ \\ \end{array} \tag{26}$$

Finally, each set of constants {*c*1, . . . , *cd*} is given by rows of *C*. Hence, each *Y<sup>i</sup>* = *w*<sup>∗</sup> · *w*<sup>T</sup> can be calculated, and the beamformers *w<sup>i</sup>* for target user signal recovery *s* T *<sup>i</sup>* <sup>=</sup> *<sup>w</sup>*<sup>T</sup> *i V*ˆ can be estimated as the singular vectors corresponding to the largest singular value in each *Y<sup>i</sup>* . The user signal space can be successfully separated from the interference space. However, sometimes, we still need to run the Gerchberg–Saxton algorithm (GSA) several

times to solve the phase invariance problem in order to make the solutions near optimal. At the *k*−th iteration, the update rule is

$$\mathbf{w}\_{i}^{\mathsf{T}}(k+1) = \begin{bmatrix} \frac{(\mathbf{w}\_{i}^{\mathsf{T}}(k)\hat{\mathcal{V}})\_{1}}{|(\mathbf{w}\_{i}^{\mathsf{T}}(k)\hat{\mathcal{V}})\_{1}|}, \dots, \frac{(\mathbf{w}\_{i}^{\mathsf{T}}(k)\hat{\mathcal{V}})\_{n}}{|(\mathbf{w}\_{i}^{\mathsf{T}}(k)\hat{\mathcal{V}})\_{n}|} \end{bmatrix}. \tag{27}$$

#### **5. Performance Evaluation**

In this section, we perform numerical simulations to test the efficiency of the analytical method for blind beamforming in a UAV-BS. We consider an uplink wireless communication scenario in which a UAV-BS serves *d* ground users in a specific urban region—for instance, a 200 m × 200 m square region. In this region, the UAV-BS hovers in the sky, while the users stay on the ground. The time-varying channel is considered according to the model in Equation (10). We assume that the UAV moves at a speed of 10 m/s [28], while the ground users maintain static. The carrier frequency is set to 28 GHz in the mmWave band. So, *fDoppler* can be easily derived. The positions of the UAV-BS and ground users are randomly distributed in the horizontal direction. For each position setting in each simulation run, we average the simulation results over 100 times.

The UAV-BS is equipped with an *m* = *m<sup>x</sup>* × *m<sup>y</sup>* element URA, while each ground user is equipped with a single antenna. In this uplink process, ground users transmit independent and randomly generated signals by using QPSK modulations with *n* snapshots. Each signal is assumed to have *L<sup>i</sup>* = 4 MPCs with the existence of LoS and the three strongest NLoS components; thus, the channel in this scenario can be regarded as a mixture of LoS and NLoS parts. We set the path-loss factors *α* = 0.95 and *β* = 2.25 in Equations (8) and (9) according to the mmWave channel measurements in Reference [37]. The small-scale factor *σ<sup>f</sup>* can be set to a small value, such as 0.01. By combining the channel model mentioned in Section 2.2 and the parameters in this section, we can simulate different channels for the tests of this method. According to detailed parameter and environment settings in this section, combined with the algorithm and models mentioned above, we can have the antenna outputs matrix *X*, signal matrix *S* and the beamforming matrix *H*. Then, we can easily evaluate the algorithm performance in this specific-simulated application scenario.

In the following simulations, additive Gaussian white noise is added to the system. We first assume that *m* = 4 × 4, and there are *δ* = 6 signals that appear in the system; four of them are transmitted by ground users (target CM signals). By exploring the number of singular near-zero values after *SVD*(*P*ˆ), we can estimate how many ground users are located in the region served by the UAV-BS. As shown in Figure 2, we tested the efficiency of the estimation of the number of ground users under different SNR and snapshot settings.

(**a**) Different SNR settings (n = 1000) (**b**) Snapshot settings (SNR = 15 dB) **Figure 2.** Estimation of the number of ground users by examining small singular values when using different SNR/snapshot settings.

Due to the existence of Gaussian white noise, the singular values corresponding to the ground users are raised above zero. However, there is still a gap between the near-zero values and other values. With the increasing SNR and number of snapshots, the gap becomes more distinct. A reasonable threshold was given in Reference [32], which the large singular values should satisfy (with a probability better than 95%)

$$
\min(\text{large singular value}) > \frac{1}{\sqrt{n}} - \frac{\delta - 0.5}{n}.\tag{28}
$$

We then tried to evaluate the performance of the analytical method in minimizing the system's MSE under different channel settings. For the LoS channel, the LoS component dominates the channel response. For the NLoS channel, all MPCs are regarded as NLoS components. The Rayleigh fading channel represents an ideal NLoS situation in which *H* ∼ *Rayleigh*(1). We used the classic MMSE method, random beamforming, and the Artificial Bee Colony (ABC) method to compare with the analytical method. As we do not know the channel *H* in a blind situation, the MMSE receiver cannot be directly employed in practice. However, here, we assume that the MMSE receiver already knows the channel. The classic MMSE will give us a weight vector

$$\mathbf{W}\_{\rm MMSE} = \mathbf{H}^{\rm H} \left( H \mathbf{H}^{\rm H} + \frac{\sigma^2}{P} \mathbf{I} \right)^{-1}, \tag{29}$$

where *σ* 2 is the system noise power, and *P* is the system signal power; *<sup>σ</sup>* 2 *P* is actually *SNR*−<sup>1</sup> on the natural scale. The random beamformer is defined as *Wrand* = *e <sup>j</sup>***Ψ**, where each element of the random phase matrix **Ψ** will give us a random phase within [0, 2*π*] [20]. The ABC method here is slightly different from the ABC beamforming process in Reference [20]; we transform it into an uplink process here. In this test, we set the food source number *Ns* = 200 with *iteration* = 500.

From Figure 3, we can observe that, for the LoS and NLoS channel settings, the analytical blind beamforming achieves better performance in minimizing the system's MSE, while the MMSE and ABC receivers do not perform as well. For the Rayleigh fading channel, the efficiency of the analytical method is close to that of the MMSE because, in this situation, noise is the dominant force affecting the signal recovery. The actual values depend on the randomness of the noise and the channel distortion. Nevertheless, the analytical blind beamforming is useful in most common scenarios.

We also tested the ACMA beamforming method under some sub-6-GHz bands. As shown in Figure 4, the results are very close to each other. The carrier frequency in our simulations mainly decides the values of the path loss and *fDoppler*, but it does not affect the performance of the method in achieving a minimal system MSE.

Figure 5 shows the constellations of the originally sent user signals and recovered user signals. It directly shows the performance of the ACMA in beamforming for signal recovery. Compared with random beamforming, the analytical method can successfully derive both the amplitude and phase information. Due to the phase ambiguity of blind beamforming methods, sometimes, there are phase differences between the original signals and the recovered ones. Operations, such as differential encoding, can be applied to cancel the phase differences. However, we mainly focused on blind beamforming for signal recovery; thus, we have omitted the details about the encoding part.

**Figure 4.** The MSE performance of the analytical method with different frequency bands.

**Figure 5.** Constellation diagrams of the received user signals when SNR = 15 dB. The upper row shows the recovered signals of four different users using the ACMA with the original constellation. The bottom row shows the signals when using random beamforming.

In Figure 6, we show the compensated channel response for *d* = 4 different users after applying analytical beamforming. We can see that the compensated channel *wH* looks the same as an impulse. This, in essence, shows that the channel effects on each transmitted user signal are nearly erased at the UAV-BS's receiver. In other words, the ACMA receivers only recover user signals (*wiHS* = *s<sup>i</sup>* ), which can be seen as directing the receiving antennas towards different users.

Then, we evaluate the system's sum rate when applying the analytical method under different settings. The system's sum rate can be described as

$$Sumrate = \sum\_{i=1}^{d} \log\_2 \left( 1 + \frac{P \cdot |\mathbf{h}\_i^{\text{H}} \mathbf{w}\_i|^2}{\sigma^2} \right), \tag{30}$$

which summarizes achievable rates of *d* ground users. Here, the signal power *P* and noise power *σ* <sup>2</sup> are still calculated on the natural scale. For better illustration, we transform them into dB in the figures. We introduced approximate beamforming as an upper bound in our tests. Approximate beamforming is regarded as an ideal beamforming in which the beam gains are zeros along the non-user directions and are significant along the user directions [20].

**Figure 6.** Compensated channel response of four different users with the ACMA in comparison with ideal beamforming when SNR = 10 dB.

In this experiment, we first evaluated the performance of these three beamforming methods with a settled noise power of *σ* <sup>2</sup> <sup>=</sup> <sup>−</sup><sup>100</sup> dBm and varying antenna numbers *<sup>m</sup>*. Then, we used a 16-element URA to evaluate the system's sum rate for varying *σ* 2 . As shown in Figure 7, under different *m* and *σ* 2 settings, the system sum rates achieved with the analytical method are close to the approximate ones; moreover, they are much better than with the random values.

As mentioned in Section 4, the Gerchberg–Saxton algorithm (GSA) is used for the phase-retrieval problem. We proceeded with some iterations of the GSA for both analytically and randomly calculated weight vectors *w<sup>i</sup>* for more accuracy. The convergence speed of GSA was tested in our simulations. As shown in Figure 8, the beamformers *w<sup>i</sup>* obtained with the analytical method can achieve a very fast convergence speed, while the ones derived from random beamforming surely need more iterations to guarantee that they converge to their optimal solutions.

(**a**) Different m (*σ* <sup>2</sup> <sup>=</sup> <sup>−</sup>100 dBm) (**b**) Different *σ* 2 (*m* = 16) **Figure 7.** The system's sum rate performance after analytical blind beamforming.

**Figure 8.** GSA convergence speed of four users using random initial points or points calculated with the ACMA when SNR = 20 dB.

After testing the performance of the analytical method in terms of the signal recovery accuracy and sum rate, we then tested its performance in terms of complexity. In Table 1, we compare it with classic non-blind methods (LS channel estimation/MMSE estimation + MMSE receiver) and the ABC method. As different methods depend on different parameters, sometimes, it is hard to directly compare the complexity in one united form. Thus, we compare these methods when they achieve comparable MSE results and introduce an equivalent complexity. In our settings, we need at least *m* > *d* and *n* > *d* <sup>2</sup> + 1, so we can transform O(*dnm*) into O(*d* 4 ) [32]. When *d* = 4, we need about *LGSA* = 20, *Ns* = 100, and *LABC* = 300 to make the CMA and ABC converge to their optimal results, so we obtain O(*LGSA*) ≈ O(*d* 2 ), O(*Ns*) ≈ O(*d* 3 ), and O(*LABC*) ≈ O(*d* 4 ). The equivalent complexity of each method will depend on the actual parameter settings.

**Table 1.** Comparison of the complexity of the methods when achieving comparable system MSE values (notations: *m*: antenna number, *d*: user number, *n*: snapshots, *Ns*: food source number, *LGSA*: iteration number of the GSA, *LABC*: iteration number of ABC).


We then tested the simulation runtimes under specific parameter settings, as shown in Figure 9. We can see that the ACMA can achieve a comparable runtime with that of the classic methods with respect to different values of *m*. However, with an increasing user number, there is a dramatic runtime increase for the ACMA, but it is still more efficient than ABC. In fact, when the user number increases, ABC needs more food sources and iterations to obtain good results, so we set *N<sup>s</sup>* = 100 and *LABC* = 300 for only the controlled variables.

(**a**) *d* = 4, *n* = 1000, *LGSA* = 20, *Ns* = 100, *LABC* = 300. (**b**) *m* = *d*, *n* = *d* <sup>2</sup> + 1, *LGSA* = 20, *Ns* = 100, *LABC* = 300. **Figure 9.** Simulation runtimes for different parameter settings.

In relation to Figure 3, there is a trade-off problem between the complexity and the system MSE. Indeed, the analytical method is more complex than the traditional MMSE method, but it achieves more accurate signal recovery beamforming when considering new features in the UAV communication channels. Accurate channel estimation is crucial for non-blind methods, but it is still difficult to acquire due to the time variance and other features [27]. From another perspective, battery endurance and energy cost are two inevitable factors in UAV-assisted networks that must still be solved, though non-blind methods usually need pilots or training sequences, which impact the bandwidth efficiency, power consumption, and battery life [38]. The computational complexity is not as difficult to deal with because chips are becoming more capable.

Through the evaluations in this section, we can observe that, when applied to a simulated UAV-BS system, the analytical blind beamforming method can achieve good MSE performance in terms of signal recovery and a reasonable system sum rate. Complex and difficult MIMO channel estimation is omitted to save costs. In other words, this method can successfully solve the MIMO UAV-BS beamforming problem without CSI. At the same time, it can achieve a reasonable computational complexity when there are not too many ground users. One possible UAV-BS application scenario can be a situation in which some specific users/target signals need higher-quality communications.

#### **6. Conclusions**

In this paper, we investigated the blind beamforming method with the ACMA at URA receivers for an MIMO UAV-BS system in mmWave bands. We firstly formulated the UAV-BS beamforming problem as a system MSE minimization problem. Then, we introduced the analytical method, ACMA, in order to derive a proper beamformer for signal recovery. Instead of relying on detailed CSI, the analytical method only explores information from antenna outputs. This method separates the user signal space from other interference spaces by using algebraic operations, such as SVD and diagonalization; then, it uses a specific structure to restrict the signal space to obtain a proper beamformer. We evaluated this method under different settings, such as different channels, different antenna settings, and different noise environments. New channel features, such as UAV mobility and time variance, were taken into consideration. The simulation results show that the analytical method exhibits good performance in terms of minimizing the system MSE, and, at the same time, it achieves a reasonable system sum rate. When the UAV-BS is not serving too many users, it achieves a complexity that is comparable with that of the MMSE method. Future research might study the overall power consumption increased by computation and transmission length. Practical experiments might be included in future work.

**Author Contributions:** Funding acquisition, project administration, supervision, K.F.; methodology, investigation, writing—original draft, P.L.; writing—review and editing, Y.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (No. 61763 018), the 03 Special Project and 5G Program of the Science and Technology Department of Jiangxi Province (No. 20193ABC03A058), the China Scholarship Council (CSC, No. 201708360150), the Key Foundation of Education Committee of Jiangxi (No. GJJ170493, No. GJJ190451), the Program of Qingjiang Excellent Young Talents in Jiangxi University of Science and Technology (JXUSTQJBJ2019004), and the cultivation project of the State Key Laboratory of Green Development and High-Value Utilization of Ionic Rare-Earth Resources in Jiangxi Province (20194AFD44003).

**Data Availability Statement:** Data are only available upon request due to restrictions regarding, e.g., privacy and ethics. The data presented in this study are available from the corresponding author upon request. The data are not publicly available due to their relation to other ongoing research.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


**Adrian Marius Deaconu <sup>1</sup> , Razvan Udroiu 2,\* and Corina- ¸Stefania Nanau <sup>1</sup>**


**Abstract:** Drones are frequently used for the delivery of materials or other goods, and to facilitate the capture and transmission of data. Moreover, drone networks have gained significant interest in a number of scenarios, such as in quarantined or isolated areas, following technical damage due to a disaster, or in non-urbanized areas without communication infrastructure. In this context, we propose a network of drones that are able to fly on a map covered by regular polygons, with a well-established mobility schedule, to carry and transfer data. Two means exist to equidistantly cover an area with points, namely, grouping the points into equilateral triangles or squares. In this study, a network of drones that fly in an aerial area divided into squares was proposed and investigated. This network was compared with the case in which the area is divided into equilateral triangles. The cost of the square drone network was lower than that of the triangular network with the same cell length, but the efficiency factors were better for the latter. Two situations related to increasing the drone autonomy using drone charging or battery changing stations were analyzed. This study proposed a Delay Tolerant Network (DTN) to optimize the transmission of data. Multiple simulation studies based on experimental flight tests were performed using the proposed algorithm versus five traditional DTN methods. A light Wi-Fi Arduino development board was used for the data transfer between drones and stations using delivery protocols. The efficiency of data transmission using single-copy and multiple-copy algorithms was analyzed. Simulation results showed a better performance of the proposed Time-Dependent Drone (TD-Drone) Dijkstra algorithm compared with the Epidemic, Spray and Wait, PRoPHET, MaxProp, and MaxDelivery routing protocols.

**Keywords:** drones; network; DTN; mobility schedule; routing algorithms; data delivery

## **1. Introduction**

Delay tolerant networks (DTNs) allow communication in environments in which frequent transmission discontinuities are present [1–4]. They have applications in numerous fields, such as space communication networks [1–3], smart cities [5], intelligent transport networks, rural networks, environmental monitoring networks, and vehicle networks [4,6]. Within DTNs, message transmission is based on the store-carry-forward paradigm [7,8]. Devices update their communication routes based on the topological changes of the network, and the mobility of the devices plays an important role [9–11].

The role of routing [12–14] in DTNs is to find the best path to send data through the network to reach the destination. The routing strategies used by DTNs are classified based on criteria such as connection type between nodes, the time at which the path for messages is established, the amount of information held by the nodes about the network, and the number of copies of a message that a node sends.

Algorithms containing different amounts of DTN information are proposed for investigation in [14]. It has been shown that the performances of these algorithms gradually increase, depending on the amount of information about the network they use. Based on

**Citation:** Deaconu, A.M.; Udroiu, R.; Nanau, C.-¸S. Algorithms for Delivery of Data by Drones in an Isolated Area Divided into Squares. *Sensors* **2021**, *21*, 5472. https://doi.org/10.3390/ s21165472

Academic Editor: Margot Deruyck

Received: 26 July 2021 Accepted: 10 August 2021 Published: 13 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the number of copies of a message, single-copy algorithms (forward based) and multiplecopy algorithms (flood based) have been investigated [15,16]. The Direct Delivery algorithm only sends the message to the destination node, and is suitable for both small- and high-mobility networks, in which the probability of meeting between the nodes is high. Comparisons between classical flood-based protocols, such as Epidemic [17], Spray and Wait [18], PRoPHET [19], and MaxProp [20], were undertaken in [15]. Maximum flow with the static approach in buffer-limited delay tolerant networks was investigated in [21,22].

In public transport networks, the connection type between nodes is based on a wellknown schedule and well-defined routes. A DTN composed of pedestrians and cyclists equipped with smart devices was investigated using routing protocols in an Opportunistic Network Environment (ONE) [23]. The MaxProp protocol yielded the best results in terms of delivery probability and average latency. DTN communications in the network storage depend on the store-carry-forward mechanism. A DTN system applied in a communication network on a railway line was found to reduce the message delivery time by 20%, depending on the schedule of trains [8]. In [24], an algorithm to search for the shortest safe path on the network with a time-dependent and edge-length danger factor was proposed. This work is suitable for the optimization of heavy trucks carrying inflammable materials, poison gas, or explosive cargo, and traveling within a city. Based on the simulation in ONE software, the Scheduling-Probabilistic Routing Protocol using History of Encounters and Transitivity (PROPHET) improves the delivery rate and optimizes the delivery delay with low overhead in DTNs for IoT applications [25].

Drone networks are now used in numerous applications across domains including topography [26], aerial observation [27], delivery [28,29], agriculture [30], communications [25,31,32], atmospheric sciences [33–35], and rescue missions [36].

Drone flight path planning can be categorized as off-line planning, on-line planning, or cooperative planning [37,38]. The relatively short flight distance of drones due to their limited battery energy [39] can be extended using drone networks.

Flight simulations using a topological map of hexagons, in which each hexagon contains one drone, have shown the potential for application in indoor rescue missions [36]. In [40], coverage path planning with UAVs was studied, addressing simple geometric flight patterns, such as back-and-forth and spiral, and more complex grid-based solutions considering full and partial information about the area of interest. The area was divided into squares and a drone was flown from a square to another. Parcel delivery missions using a drone were simulated based on heuristic flight path planning (HFPP) and other routing algorithms [28]. The results showed that HFPP delivers up to 33% more data packets compared with Encounter-Based Routing and Epidemic routing protocols.

Vehicular ad hoc networks (VANETs), which are a type of mobile ad hoc network (MANET), can be used in an intelligent transport system. VANETs allow the mobile vehicles to establish three main categories of communication: vehicle to vehicle, vehicle to infrastructure, and infrastructure to infrastructure. A specific application of VANET applied to a drone network allows messages to be sent via wireless links [32]. A bio-inspired coordination protocol for a drone flying ad-hoc network (FANET) used for agriculture applications [30] has been used in an ad hoc simulator for a preliminary analysis of the feasibility of drone network design. The performance of the branch and bound searchbased mode selection (BBS-MS) for drone-based air-to-ground wireless networks was investigated in [31].

In this study, we investigated six routing algorithms in a network of drones with a mobility schedule to ensure communication between isolated areas or in areas with technical damage. We then qualitatively compared the algorithms. The major contributions of this study are as follows:


The remainder of the paper is organized as follows. In Section 2, a drone network architecture based on the mission profile, and drone sensing and communication, are proposed. Then, the algorithms and protocols for the delivery of data using drones are proposed, described, and analyzed. In Section 3, simulation results are summarized and discussed. Finally, conclusions and possible future work are presented in Section 4.

#### **2. Materials and Methods**

#### *2.1. Drone Network Architecture and Communication*

A network map is essential for both drone flight control and simultaneous localization and mission tasks. As noted in [36], three forms of map generation exist: metric, topological, and hybrid maps. A metric map is represented as a grid, geometric, or feature map. Topological maps are represented by graphs comprised of nodes and edges, where nodes represent places, and the edges represent the paths between the nodes [36]. A hybrid map consists of small metric map locations in nodes. These nodes are connected by edges, which are the paths between the metric maps [36]. In this work, a topological map was proposed.

The 2D network surface is intended to be covered with equidistantly distributed points. These points form regular polygons covering the 2D surface. Below, we prove that the only possible means of covering the 2D surface in this manner is by using equilateral triangles or squares. First, we prove the following lemma:

**Lemma 1.** *There are 3 ways to cover a 2D surface using regular polygons: equilateral triangles, squares, or regular hexagons.*

**Proof of Lemma 1.** The sum of the degrees of the angles of a polygon with *n* vertices (*n* ≥ 3) is 180◦ · (*n* − 2). Each angle of a regular polygon has [360◦ · (*n* − 2)]/*n* degrees. If a point of the surface is a vertex of a polygon, then it is a vertex for *m* polygons (*m* ≥ 3) around this point (Figure 1). Thus, at such a point we have:

$$m \cdot \frac{180^{\circ} \cdot (n - 2)}{n} = 360^{\circ} \Leftrightarrow \frac{2 \cdot n}{n - 2} = m \tag{1}$$

It results that 2 · *n*/(*n* − 2) = 2 + 4/(*n* − 2) is an integer value greater or equal to 3. This implies that *n* − 2 is a divisor of 4 and, because *n* ≥ 3, it follows that *n* − 2 is one of the values 1, 2, or 4, which is equivalent to the fact that *n* is 3, 4, or 6, and m is 6, 4 or, respectively, 3. This means that the 2D surface can only be covered with equilateral triangles, squares, or regular hexagons.

**Theorem 1.** *There are two ways to equidistantly cover a 2D surface with points.*

**Proof of Theorem 1.** Using the result from Lemma 1 and the fact that the points are considered equidistant on the 2D surface, it follows that there are only two ways to cover

the 2D surface: with equilateral triangles or with squares, because using regular hexagons is impossible (the vertices of the hexagons are not equidistantly positioned). *Sensors* **2021**, *21*, x FOR PEER REVIEW 4 of 21

**Figure 1.** Regular polygons sharing the same vertex. **Figure 1.** Regular polygons sharing the same vertex.

It results that 2 ⋅ *n* /(*n* − 2) = 2 + 4/(*n* − 2) is an integer value greater or equal to 3. This implies that *n* − 2 is a divisor of 4 and, because *n* ≥ 3, it follows that *n* − 2 is one of the values 1, 2, or 4, which is equivalent to the fact that *n* is 3, 4, or 6, and m is 6, 4 or, respectively, 3. This means that the 2D surface can only be covered with equilateral triangles, squares, or regular hexagons. **Theorem 1**. *There are two ways to equidistantly cover a 2D surface with points.*  **Proof of Theorem 1.** Using the result from Lemma 1 and the fact that the points are considered equidistant on the 2D surface, it follows that there are only two ways to cover A network to cover a surface with hexagon cells containing three equilateral triangles was proposed in [41]. In the current paper, a network is proposed to cover a similar surface with square cells containing two operational squares (Figure 2). Each operational cell is covered by a drone. A battery charging/changing dock, which is shared by two drones, is placed in the center of each cell. The docks are able to automatically change/recharge drones, without manual intervention, allowing fully autonomous drone management. The mission profile of each drone consists of several phases or steps: engines start, take off, climb at the cruise altitude, cruise, hovering and data exchange, descent, landing, and engines shut-down. The cruise phase of the flight mission profile consists of four segments (Figure 3). Charging or changing the drone battery is a necessary step at the end of the flight mission. *Sensors* **2021**, *21*, x FOR PEER REVIEW 5 of 21

**Figure 2.** Square cell with two drones. **Figure 2.** Square cell with two drones.

change points via wireless links.

**Table 1.** Drone performance [43].

**Figure 3.** Square shape flight mission profile of the drone.

A drone with a quadcopter configuration [42], i.e., a DJI Mavic 2 Pro (DJI, Shenzhen, China) with a size of 214 × 91 × 84 mm (length × width × height), and takeoff weight of 905 g, was considered in this study. The performance characteristics of this drone are presented in Table 1 [43]. All of the drones of the proposed squares network operate in a pre-programmed manner. The drone network communicates with data ex-

The endurance of drones can be improved using various methods, such as changing the battery [44,45], charging the battery via wires [46,47], wireless recharging [48], solar cells [49,50], laser-beam in-flight recharging [49,51], and tethered drones [49]. An average time of 90 min is needed to fully charge an empty battery. An automated means of charging a battery via a wire can be performed using a charging platform installed on the ground and a drone retrofit-kit mounted on the drone [46,47]. Thus, the landing gear of the drone is connected by touch with the charging platform after the drone lands, and

**Parameter Value** 

Max flight time (no wind) 31 min (at a consistent 25 km/h) Max flight distance (no wind) 18 km (at a consistent 50 km/h) Drone battery 3850 mAh, 1800 mA, 3.83 V

Max ascent/descent speed 4 m/s; 3 m/s

**Figure 3.** Square shape flight mission profile of the drone. **Figure 3.** Square shape flight mission profile of the drone.

**Figure 2.** Square cell with two drones.

A drone with a quadcopter configuration [42], i.e., a DJI Mavic 2 Pro (DJI, Shenzhen, China) with a size of 214 × 91 × 84 mm (length × width × height), and takeoff weight of 905 g, was considered in this study. The performance characteristics of this drone are presented in Table 1 [43]. All of the drones of the proposed squares network operate in a pre-programmed manner. The drone network communicates with data ex-A drone with a quadcopter configuration [42], i.e., a DJI Mavic 2 Pro (DJI, Shenzhen, China) with a size of 214 × 91 × 84 mm (length × width × height), and takeoff weight of 905 g, was considered in this study. The performance characteristics of this drone are presented in Table 1 [43]. All of the drones of the proposed squares network operate in a pre-programmed manner. The drone network communicates with data exchange points via wireless links.

change points via wireless links. **Table 1.** Drone performance [43].


Max flight distance (no wind) 18 km (at a consistent 50 km/h) Drone battery 3850 mAh, 1800 mA, 3.83 V The endurance of drones can be improved using various methods, such as changing the battery [44,45], charging the battery via wires [46,47], wireless recharging [48], solar cells [49,50], laser-beam in-flight recharging [49,51], and tethered drones [49]. An average time of 90 min is needed to fully charge an empty battery. An automated means of charging a battery via a wire can be performed using a charging platform installed on the ground and a drone retrofit-kit mounted on the drone [46,47]. Thus, the landing gear The endurance of drones can be improved using various methods, such as changing the battery [44,45], charging the battery via wires [46,47], wireless recharging [48], solar cells [49,50], laser-beam in-flight recharging [49,51], and tethered drones [49]. An average time of 90 min is needed to fully charge an empty battery. An automated means of charging a battery via a wire can be performed using a charging platform installed on the ground and a drone retrofit-kit mounted on the drone [46,47]. Thus, the landing gear of the drone is connected by touch with the charging platform after the drone lands, and charging starts automatically. The main disadvantage of this system is that the drone is locked on the ground during the charging of the battery.

of the drone is connected by touch with the charging platform after the drone lands, and An automatic battery changing and recharging system was investigated in [44]. The battery is automatically recharged after it is changed at the station. The electricity needed at each station is able to be provided by a solar panel that charges a battery located at the station. The changing time depends on the efficiency of the changing mechanism, and varies between 15 s [44] and 60 s [45]. In the current study, a maximum changing time before the drone is ready to take off of 60 s was considered.

Four flight tests were performed based on the square-shaped flight mission in the Brasov area of Romania. The flight tests were performed at a temperature of 2 ◦C, humidity of 69%, and wind speed of 3.5 km/h. A Samsung S9 smartphone device, on which DJI Go 4 app software (DJI, Shenzhen, China) was installed, connected to a DJI remote controller, was used to program the flight mission segments and to remotely control the drone (Figure 4). A flight altitude of 30 m was chosen [52] and the cruise flight distance of each drone was 12,000 m.

locked on the ground during the charging of the battery.

time before the drone is ready to take off of 60 s was considered.

**Figure 4.** Square-shaped flight mission profile. **Figure 4.** Square-shaped flight mission profile.

of each drone was 12,000 m.

The theoretical flight times calculated based on the drone specifications (max. ascent, descent speeds, and cruise speed) are not realistic for simulations. The main factors that influence the experimental flight times are acceleration and deceleration of the drone, and wind speed. The mean flight time for each flight segment was calculated based on the flight tests (Table 2). Moreover, the average calculated cruise speed of the drone obtained from the flight tests was 12.93 m/s. The theoretical flight times calculated based on the drone specifications (max. ascent, descent speeds, and cruise speed) are not realistic for simulations. The main factors that influence the experimental flight times are acceleration and deceleration of the drone, and wind speed. The mean flight time for each flight segment was calculated based on the flight tests (Table 2). Moreover, the average calculated cruise speed of the drone obtained from the flight tests was 12.93 m/s.

charging starts automatically. The main disadvantage of this system is that the drone is

An automatic battery changing and recharging system was investigated in [44]. The battery is automatically recharged after it is changed at the station. The electricity needed at each station is able to be provided by a solar panel that charges a battery located at the station. The changing time depends on the efficiency of the changing mechanism, and varies between 15 s [44] and 60 s [45]. In the current study, a maximum changing

Four flight tests were performed based on the square-shaped flight mission in the Brasov area of Romania. The flight tests were performed at a temperature of 2 °C, humidity of 69%, and wind speed of 3.5 km/h. A Samsung S9 smartphone device, on which DJI Go 4 app software (DJI, Shenzhen, China) was installed, connected to a DJI remote controller, was used to program the flight mission segments and to remotely control the drone (Figure 4). A flight altitude of 30 m was chosen [52] and the cruise flight distance

**Table 2.** Flight test results. **Table 2.** Flight test results.


The percentage of the remaining drone battery obtained at the end of the flight mission was 16% for the rectangular cell. The charging time for the drone battery was 79 min. A safety multiplier of 1.1 was applied to obtain the considered charging time, and the resulting time was used in the simulations. Thus, 87 charging minutes were considered for the rectangular cell flight. The average values of the flight times for each segment were used as input parameters for the simulation of the drone networks in the The percentage of the remaining drone battery obtained at the end of the flight mission was 16% for the rectangular cell. The charging time for the drone battery was 79 min. A safety multiplier of 1.1 was applied to obtain the considered charging time, and the resulting time was used in the simulations. Thus, 87 charging minutes were considered for the rectangular cell flight. The average values of the flight times for each segment were used as input parameters for the simulation of the drone networks in the DTN algorithms.

DTN algorithms. Eight high-resolution and two infrared sensors were used on the DJI Mavic 2 Pro. These sensors enabled omnidirectional obstacle sensing, to determine the relative speed and distance between the drone and the object, and to ensure good stability in forward Eight high-resolution and two infrared sensors were used on the DJI Mavic 2 Pro. These sensors enabled omnidirectional obstacle sensing, to determine the relative speed and distance between the drone and the object, and to ensure good stability in forward and hovering flight. Left, right, up, down, forward, and backward obstacle sensing were used.

A NodeMCU Lua Wi-Fi, V3, ESP-12E, CP2102 Wi-Fi Arduino development board (Espressif Systems, Shanghai, China) was used for the data transfer between the drones and the stations. The main characteristics of the NodeMCU are shown in Table 3. The data were stored on a micro-SD card. The Wi-Fi board and micro-SD card module are very lightweight and have very low power consumption. The experimental layout, consisting of the NodeMCU Lua Wi-Fi board and the micro-SD card module mounted on a breadboard, is shown in Figure 5a. Arduino code was used to program the Wi-Fi boards.

**Table 3.** Wi-Fi NodeMCU main characteristics [53].


used.

mounted on it.

**Table 3.** Wi-Fi NodeMCU main characteristics [53].

**Parameter Value** 

Dimensions (L × W) 48 mm × 25 mm Operating temperature −40 °C to + 125 °C Weight 8 g

ESP8266 chip 26 MHz, 4 MB flash, 160 KB RAM

A protective mounting case (Figure 5b) for the Wi-Fi boards was designed using the SolidWorks version 2016 software (Dassault Systèmes, Waltham, MA, USA), and then 3D printed using material extrusion technology from PLA (BCN3D Technologies, Barcelona, Spain) on a BCN3D Sigma R19 printer (BCN3D Technologies, Barcelona, Spain).

and hovering flight. Left, right, up, down, forward, and backward obstacle sensing were

A NodeMCU Lua Wi-Fi, V3, ESP-12E, CP2102 Wi-Fi Arduino development board (Espressif Systems, Shanghai, China) was used for the data transfer between the drones and the stations. The main characteristics of the NodeMCU are shown in Table 3. The data were stored on a micro-SD card. The Wi-Fi board and micro-SD card module are very lightweight and have very low power consumption. The experimental layout, consisting of the NodeMCU Lua Wi-Fi board and the micro-SD card module mounted on a breadboard, is shown in Figure 5a. Arduino code was used to program the Wi-Fi boards.

The Wi-Fi board, micro-SD module, SD card, connection wires, and power cable weighed a total of 21.8 g, thus representing an increase of 2.4% in the total weight of the drone. A protective mounting case (Figure 5b) for the Wi-Fi boards was designed using the SolidWorks version 2016 software (Dassault Systèmes, Waltham, MA, USA), and then 3D printed using material extrusion technology from PLA (BCN3D Technologies, Barcelona, Spain) on a BCN3D Sigma R19 printer (BCN3D Technologies, Barcelona, Spain).

*2.2. Algorithms and Protocols for Delivery of Data Using Drones*  The DTN was modeled with a graph having fixed and mobile nodes. There are no The Wi-Fi board, micro-SD module, SD card, connection wires, and power cable weighed a total of 21.8 g, thus representing an increase of 2.4% in the total weight of the drone.

#### connections between the fixed nodes, so there is no possibility of direct data transmission because the considered distance between the nearest two nodes is 3000 m. Network *2.2. Algorithms and Protocols for Delivery of Data Using Drones*

connections are provided by mobile nodes (drones), but their condition is not always the same; they have periods when they are active and periods when they are inactive. This means that there is not always an end-to-end path available between any two nodes in the graph. The DTN was modeled with a graph having fixed and mobile nodes. There are no connections between the fixed nodes, so there is no possibility of direct data transmission because the considered distance between the nearest two nodes is 3000 m. Network connections are provided by mobile nodes (drones), but their condition is not always the same; they have periods when they are active and periods when they are inactive. This means that there is not always an end-to-end path available between any two nodes in the graph.

We tested five well-known routing algorithms for DTNs (Epidemic, Spray and Wait, PRoPHET, MaxProp, and MaxDelivery [54]), and a newly proposed TD-Drone Dijkstra approach, on the square-shaped network shown in Figure 6. We performed tests on the network by choosing random sources and random destinations that could be located on any node (marked as a gray circle or a gray rectangle). We considered the drones worked each day from 7:00 a.m. to 6:00 p.m. and the messages may leave a source node between 7:00 a.m. and 5:00 p.m. We considered 1000 messages randomly sent within this interval of time.

Epidemic is the basic form of a flood-based routing protocol: when two nodes meet, they identify the packages that the other node has and transfer the packages that it does not have. At the end of the process, the two nodes have the same content in the buffer. This process is repeated each time two nodes come into contact. When a node has a copy of a message, it waits to meet the destination. In this case, the resource consumption is high but, in a high mobility network, the delay of message transmission is small. In the current network configuration, the algorithm produces poor results due to the low number of contacts between nodes.

Spray and Wait is an algorithm with two phases, one for sending messages (spray) and one to wait for the contact with the destination node (wait). This algorithm circulates in two variants—standard and binary—depending on the number of spread copies of the message. It acts in the same manner as Epidemic, with an important difference: the number of spread copies is constant. The spray phase of the standard approach consists of spraying

L copies of the message by the source node itself. The spray phase of the binary approach consists of spraying half of the number of copies to a meeting node. In this case, not only does the source spray messages, but also every node that has more than one copy. The nodes that have only one copy enter the wait phase. This algorithm has the disadvantage that nodes must keep track of other nodes' movement, but the advantage is that the level of flooding is limited. not have. At the end of the process, the two nodes have the same content in the buffer. This process is repeated each time two nodes come into contact. When a node has a copy of a message, it waits to meet the destination. In this case, the resource consumption is high but, in a high mobility network, the delay of message transmission is small. In the current network configuration, the algorithm produces poor results due to the low number of contacts between nodes.

Epidemic is the basic form of a flood-based routing protocol: when two nodes meet, they identify the packages that the other node has and transfer the packages that it does

We tested five well-known routing algorithms for DTNs (Epidemic, Spray and Wait, PRoPHET, MaxProp, and MaxDelivery [54]), and a newly proposed TD-Drone Dijkstra approach, on the square-shaped network shown in Figure 6. We performed tests on the network by choosing random sources and random destinations that could be located on any node (marked as a gray circle or a gray rectangle). We considered the drones worked each day from 7:00 a.m. to 6:00 p.m. and the messages may leave a source node between 7:00 a.m. and 5:00 p.m. We considered 1000 messages randomly

*Sensors* **2021**, *21*, x FOR PEER REVIEW 8 of 21

sent within this interval of time.

**Figure 6.** Squares network considered for experiments. **Figure 6.** Squares network considered for experiments.

Spray and Wait is an algorithm with two phases, one for sending messages (spray) and one to wait for the contact with the destination node (wait). This algorithm circulates in two variants—standard and binary—depending on the number of spread copies of the message. It acts in the same manner as Epidemic, with an important difference: the number of spread copies is constant. The spray phase of the standard approach consists of spraying L copies of the message by the source node itself. The spray phase of the binary approach consists of spraying half of the number of copies to a meeting node. In this case, not only does the source spray messages, but also every node that has more PRoPHET is similar to Epidemic, with the exception that it uses information from the buffer of the other node to update its predictability vector. Each node calculates the predictability of the message delivery and sends the message onward only if the contact node has higher predictability than its own. The disadvantage of this approach is the relationship between the overhead ratio and the number of nodes—as the number of nodes increases, the overhead ratio increases [15]. This protocol is known for the complexity of its forwarding strategy. Thus, it consumes a significant quantity of resources to process and store historical values. This approach is feasible for networks with high computation and infrastructure capabilities.

than one copy. The nodes that have only one copy enter the wait phase. This algorithm has the disadvantage that nodes must keep track of other nodes' movement, but the advantage is that the level of flooding is limited. PRoPHET is similar to Epidemic, with the exception that it uses information from the buffer of the other node to update its predictability vector. Each node calculates the predictability of the message delivery and sends the message onward only if the contact node has higher predictability than its own. The disadvantage of this approach is the relationship between the overhead ratio and the number of nodes—as the number of MaxProp is an algorithm based on prioritizing packet transmission and discarding. The packets in the queue are divided into two categories: those below the "n" hop threshold (up to that point), and those above this threshold. Newer packages that have not traveled too much are considered a priority and the guarantee that they will reach their destination is considered to be high. This algorithm also requires high computation and infrastructure capabilities. This protocol has low performance when nodes have small buffer sizes because of the adaptive threshold calculation, but gives better performance with a larger buffer size. It has a so-called slow start problem, because, in the case of a big network, it may take a very long time before each node receives the delivery predictability of other nodes because of the disconnecting nature of the networks, as shown in [20].

MaxDelivery [54] is an algorithm based on prioritizing message delivery using an appropriate buffer management strategy that consists of a forwarding, dropping, and buffer-cleaning mechanism.

Next, we propose our time-dependent Dijkstra algorithm to find a route between two nodes in our drone network. Dijkstra's algorithm is used to find the shortest path connecting two nodes in a network [55]. Our routing problem can be modeled using a time-dependent oriented network [56] defined as follows.

**Definition 1.** *A triple G = (V, A, f) is called a time-dependent oriented network, where V is a set of vertices, A* ⊆ *V* × *V is a set of arcs, and f is the time dependency function defined on each arc, f: A* × *T* → *T, T* ⊆ *R+ is called the time set and when moving on the arc a = (u, v)* ∈ *A from node u to node v, f(a, t)* ∈ *T is the moment of arrival at node v if node u is left at the moment t* ∈ *T. Of course, f(a, t) > t, for each arc a* ∈ *A and for any moment t* ∈ *T.*

In our problem, V is the equidistantly distributed set of 2D points, A is the set of connections between these points ensured by drones, and T = {0, 1, 2, . . . } is the set of counting seconds in a day. The value of f(a, t) must be computed as quickly as possible. Thus, for each arc a = (u, v), the arrival moments of drones at node v are maintained in order to enable a binary search to be undertaken when f(a, t) is calculated at arc a and time t. The arrival moments of drones for all arcs are pre-calculated (once before starting to use the algorithm) because an exact schedule of drones is known based on each drone's starting second in a day, its travel on each arc, its data transfer at nodes, and its wireless charge/battery change time.

For a given source node s ∈ V, the function dists: V × T → N is introduced, where dists(u, ts) is the distance in seconds from s to u starting from the moment ts, i.e., it is the minimum number of seconds needed to go from node s to node u ∈ V on the arcs of A if source s is left at moment ts. Of course, dists(s, ts) = 0 for every ts ∈ T. If node u is not reachable (accessible) from s, then dists(u, ts) = +∞. For a given destination node d, our problem is to determine dists(d, ts) at a given time ts. The pseudo-code for the time-dependent Dijkstra algorithm is presented in Figure 7. *Sensors* **2021**, *21*, x FOR PEER REVIEW 10 of 21

**Figure 7.** Pseudo-code of TD-Drone Dijkstra algorithm. **Figure 7.** Pseudo-code of TD-Drone Dijkstra algorithm.

(from d to s), the route must be reversed at the end.

Q is a priority queue. This means that, at any moment, the nodes from Q are sorted

are used to determine the route from s to d. If a node v was reached, then this was done using the arc (p(v), v), where node p(v) is called the predecessor of v. After the timedependent Dijkstra algorithm is executed using the staring moment ts, the route from s to d (if it exists, i.e., dists(d) < +∞) is determined using the algorithm presented in Figure 8. Because the nodes of the route are found by the above algorithm in inverse order

Q is a priority queue. This means that, at any moment, the nodes from Q are sorted in ascending order according to their distance in seconds from s. At the end of the algorithm, if destination d was reached, the values kept in the so-called predecessor vector p are used to determine the route from s to d. If a node v was reached, then this was done using the arc (p(v), v), where node p(v) is called the predecessor of v. After the time-dependent Dijkstra algorithm is executed using the staring moment ts, the route from s to d (if it exists, i.e., dists(d) < +∞) is determined using the algorithm presented in Figure 8. Because the nodes of the route are found by the above algorithm in inverse order (from d to s), the route must be reversed at the end. *Sensors* **2021**, *21*, x FOR PEER REVIEW 11 of 21


**Figure 8.** Pseudo-code for the construction of the route after the execution of the TD-Drone Dijkstra algorithm. **Figure 8.** Pseudo-code for the construction of the route after the execution of the TD-Drone Dijkstra algorithm.

For each data file that has to be delivered, a "json" file is attached that stores all of the information needed to transfer the file from the source to the destination: delivery type (Dijkstra, Epidemic etc.); file information (name, size, time-to-leave); route information (node codes: stations and drones); etc. The json files that store message information for the algorithms Epidemic, Spray and Wait, PRoPHET, MaxProp, and MaxDelivery implemented in the ONE environment are similar to those considered in [41]. Each route starts and ends with a station id. Each station id (except for the destination) is followed by a drone id, and each drone id is followed by a station id. When a drone arrives at a station, a transfer is initiated between the drone and the station. The drone transfers to the station all of the files that have the station's id in the attached json files. After the drone transfer to the station is completed, the station transfers to the drone all of the files that have the drone's id in the json file. For each data file that has to be delivered, a "json" file is attached that stores all of the information needed to transfer the file from the source to the destination: delivery type (Dijkstra, Epidemic etc.); file information (name, size, time-to-leave); route information (node codes: stations and drones); etc. The json files that store message information for the algorithms Epidemic, Spray and Wait, PRoPHET, MaxProp, and MaxDelivery implemented in the ONE environment are similar to those considered in [41]. Each route starts and ends with a station id. Each station id (except for the destination) is followed by a drone id, and each drone id is followed by a station id. When a drone arrives at a station, a transfer is initiated between the drone and the station. The drone transfers to the station all of the files that have the station's id in the attached json files. After the drone transfer to the station is completed, the station transfers to the drone all of the files that have the drone's id in the json file.

The range of the Wi-Fi boards was tested. The connections between the boards and file transfer were performed at a distance of up to 85 m with no obstacles in between. In Figure 9, the console output of the Wi-Fi board is presented. The following steps are executed: setup (upload speed is set to 115,200, MAC address is obtained), the connection between two Wi-Fi boards is established, the file transfer is performed, and, finally, the Wi-Fi boards are disconnected. The range of the Wi-Fi boards was tested. The connections between the boards and file transfer were performed at a distance of up to 85 m with no obstacles in between. In Figure 9, the console output of the Wi-Fi board is presented. The following steps are executed: setup (upload speed is set to 115,200, MAC address is obtained), the connection between two Wi-Fi boards is established, the file transfer is performed, and, finally, the Wi-Fi boards are disconnected.

In our model, when transferring files, the drone hovers over the station at a height of 30 m, which is significantly less than the maximum distance obtained in range tests. The transfer speed, including writing on and reading from the SD card, was also tested at a distance of 30 m. An average of 5.81 Mbps was obtained, which means that a file of 10 MB was transferred in 13.77 s.


**Figure 9.** The console output of the Wi-Fi board. **Figure 9.** The console output of the Wi-Fi board.

#### In our model, when transferring files, the drone hovers over the station at a height **3. Simulation Results and Discussion**

of 30 m, which is significantly less than the maximum distance obtained in range tests. The transfer speed, including writing on and reading from the SD card, was also tested at a distance of 30 m. An average of 5.81 Mbps was obtained, which means that a file of 10 MB was transferred in 13.77 s. **3. Simulation Results and Discussion**  The proposed network of drones can be applied in various scenarios, such as in remote quarantined or isolated areas, following technical damage due to a disaster (e.g., an earthquake), or in non-urbanized areas without electricity access or communication infrastructure. For instance, a remote quarantined zone (e.g., due to the COVID-19 pan-The proposed network of drones can be applied in various scenarios, such as in remote quarantined or isolated areas, following technical damage due to a disaster (e.g., an earthquake), or in non-urbanized areas without electricity access or communication infrastructure. For instance, a remote quarantined zone (e.g., due to the COVID-19 pandemic), in which buildings are separated by a safe distance, is considered. In this scenario, each building may be a house in which patients are isolated, a warehouse storing food or drugs, a laboratory in which medical tests are performed for patients, or a location at which doctors are working (isolated from patients). Data packages (medical images, tests, results, prescriptions from doctors, etc.) must be sent between these buildings. The communication between these buildings can be achieved by drones organized in a square-shaped network.

demic), in which buildings are separated by a safe distance, is considered. In this scenario, each building may be a house in which patients are isolated, a warehouse storing food or drugs, a laboratory in which medical tests are performed for patients, or a location at which doctors are working (isolated from patients). Data packages (medical im-To validate the proposed method, simulations were performed using the drone network maps as a collection of squares. Most of the simulations were performed using the Java-based simulator ONE and its provided facilities, using an ASUS ROG GL752VW-T4015D laptop with Intel® Core™ i7-6700HQ 2.60GHz processor (Asus, Taipei City, Tai-

ages, tests, results, prescriptions from doctors, etc.) must be sent between these build-

wan), and 8 GB of RAM. The routing protocols used for the simulation within the ONE simulator were Epidemic, Spray and Wait, PRoPHET, MaxProp, and MaxDelivery.

The following main steps were used for implementation of the ONE scenarios. The first step consists of defining the map (Figure 10) in wkt file format, in which the coordinates of all the points on the map, including points that establish the route of each drone, are defined. The initializations of the algorithms that define the mobility of drones consist of establishing the initial positions of drones and the recharging/changing points, associating each drone with a recharger/changing point, establishing the stationary points for data transfer, and defining the route of each drone. The final step is the establishment of the simulation parameters, as shown in Table 4. The time parameters, such as the travel autonomy time, the hovering time for the transfer points, and the parking time at the charging or changing points, were established based on the experimental flight tests of the DJI Mavic 2 Pro drone.

**Table 4.** Simulation parameters for the square-shaped flight mission.


The proposed time-dependent Dijkstra variant was implemented in Visual C++

The delivery rate and latency metrics were used to measure the performance of all six routing protocols analyzed in this paper. The delivery rate is determined as a ratio between the number of delivered messages and the number of created messages. The latency is the average time needed for a message to reach the destination starting from the

Detailed results of the comparison between changing the battery and charging the battery cases are presented in Table 5. The results obtained in this paper by simulation

In the case of the squares network and battery changing, the values of the delivery

**Battery Changing** 

**Battery Charging** 

**Delivery Rate Latency** (**hours**)

**Squares Triangular Squares Triangular Squares Triangular Squares Triangular** 

2017 programming language. The application has about 1100 lines of C++ source code.

tions: squares with battery charging, and squares with battery changing. For each case, the same 1000 route simulations considered in the ONE experiments were executed.

**Figure 10.** Design of the squares network in the ONE simulator. rate were within the range of 0.166 to 0.646 for the routing protocols Epidemic, Spray **Figure 10.** Design of the squares network in the ONE simulator.

**Table 5.** Efficiency factors in drone network.

Epidemic 0.166 0.209 0.135 0.146 0.81 0.72 2.28 2.13 Spray and Wait 0.211 0.179 0.141 0.156 0.52 0.56 1.75 1.92 PRoPHET 0.594 0.762 0.143 0.319 0.61 0.52 2.28 2.49 MaxProp 0.646 0.743 0.135 0.261 0.52 0.47 1.72 1.90 MaxDelivery 0.203 0.271 0.139 0.160 1.08 0.71 2.11 1.80 TD-Drone Dijkstra 0.954 0.973 0.540 0.664 0.43 0.45 1.69 1.48

**Battery Charging** 

source (departure node).

**Battery Changing** 

**Algorithm** 

were compared with those obtained in [41].

The proposed time-dependent Dijkstra variant was implemented in Visual C++ 2017 programming language. The application has about 1100 lines of C++ source code.

The application written in C++ was executed for each of the two considered situations: squares with battery charging, and squares with battery changing. For each case, the same 1000 route simulations considered in the ONE experiments were executed.

The delivery rate and latency metrics were used to measure the performance of all six routing protocols analyzed in this paper. The delivery rate is determined as a ratio between the number of delivered messages and the number of created messages. The latency is the average time needed for a message to reach the destination starting from the source (departure node).

Detailed results of the comparison between changing the battery and charging the battery cases are presented in Table 5. The results obtained in this paper by simulation were compared with those obtained in [41].


*Sensors* **2021**, *21*, x FOR PEER REVIEW 15 of 21

**Table 5.** Efficiency factors in drone network.

In the case of the squares network and battery changing, the values of the delivery rate were within the range of 0.166 to 0.646 for the routing protocols Epidemic, Spray and Wait, PRoPHET, MaxProp, and MaxDelivery. The best delivery rate was 0.954 for the TD-Drone Dijkstra protocol (Figure 11). The worst result in terms of latency was obtained for MaxDelivery, and TD-Drone Dijkstra's average latency was the best, as expected (Figure 12). and Wait, PRoPHET, MaxProp, and MaxDelivery. The best delivery rate was 0.954 for the TD-Drone Dijkstra protocol (Figure 11). The worst result in terms of latency was obtained for MaxDelivery, and TD-Drone Dijkstra's average latency was the best, as expected (Figure 12).

The delivery rate in the case of drone battery charging was between 0.135 and 0.143 for the routing protocols Epidemic, Spray and Wait, PRoPHET, MaxProp, and MaxDelivery. The maximum delivery rate was 0.54 and was obtained for the TD-Drone Dijkstra protocol. These low values were obtained because of the large delay due to the battery charge. The worst results for latency were obtained for PRoPHET and Epidemic, and TD-Drone Dijkstra's latency was also the best. The delivery rate in the case of drone battery charging was between 0.135 and 0.143 for the routing protocols Epidemic, Spray and Wait, PRoPHET, MaxProp, and MaxDelivery. The maximum delivery rate was 0.54 and was obtained for the TD-Drone Dijkstra protocol. These low values were obtained because of the large delay due to the battery charge. The worst results for latency were obtained for PRoPHET and Epidemic, and TD-Drone Dijkstra's latency was also the best.

**Figure 11.** Delivery rate in the squares drone network. **Figure 11.** Delivery rate in the squares drone network.

Battery changing Battery charging

and 64 routes in the case of battery recharge).

0 0.5 1 1.5 2 2.5

Average

latency [hours]

**Figure 12.** Average latency in the squares drone network.

The comparison between average latencies was performed on the routes where the

The best delivery rate in all cases was given by the TD-Drone Dijkstra algorithm because its buffer load was the lowest (the data was loaded only in the stations and drones belonging to the calculated route). The latency in the case of Dijkstra's algorithm was high because it can deliver most of the packages. Generally, the other algorithms can de-

The comparison between average latencies was performed on the routes where the delivery was successful for all of the algorithms (66 routes for the case of battery change and 64 routes in the case of battery recharge). The comparison between average latencies was performed on the routes where the delivery was successful for all of the algorithms (66 routes for the case of battery change and 64 routes in the case of battery recharge).

*Sensors* **2021**, *21*, x FOR PEER REVIEW 15 of 21

and TD-Drone Dijkstra's latency was also the best.

Battery changing Battery charging

**Figure 11.** Delivery rate in the squares drone network.

pected (Figure 12).

0 0.2 0.4 0.6 0.8 1

Delivery rate

and Wait, PRoPHET, MaxProp, and MaxDelivery. The best delivery rate was 0.954 for the TD-Drone Dijkstra protocol (Figure 11). The worst result in terms of latency was obtained for MaxDelivery, and TD-Drone Dijkstra's average latency was the best, as ex-

The delivery rate in the case of drone battery charging was between 0.135 and 0.143 for the routing protocols Epidemic, Spray and Wait, PRoPHET, MaxProp, and MaxDelivery. The maximum delivery rate was 0.54 and was obtained for the TD-Drone Dijkstra protocol. These low values were obtained because of the large delay due to the battery charge. The worst results for latency were obtained for PRoPHET and Epidemic,

**Figure 12.** Average latency in the squares drone network.

**Figure 12.** Average latency in the squares drone network.

The best delivery rate in all cases was given by the TD-Drone Dijkstra algorithm because its buffer load was the lowest (the data was loaded only in the stations and drones belonging to the calculated route). The latency in the case of Dijkstra's algorithm was high because it can deliver most of the packages. Generally, the other algorithms can de-The best delivery rate in all cases was given by the TD-Drone Dijkstra algorithm because its buffer load was the lowest (the data was loaded only in the stations and drones belonging to the calculated route). The latency in the case of Dijkstra's algorithm was high because it can deliver most of the packages. Generally, the other algorithms can deliver only on the shorter distances due to buffer restrictions. It is known that Dijkstra's algorithm results in the shortest path, and, therefore, the best latency. *Sensors* **2021**, *21*, x FOR PEER REVIEW 16 of 21 liver only on the shorter distances due to buffer restrictions. It is known that Dijkstra's algorithm results in the shortest path, and, therefore, the best latency.

> The delivery rate was better for the triangles drone network (Figures 13 and 14), with the exception of the Spray and Wait algorithm in the battery changing case. For PRoPHET, Max Prop, and TD-Drone Dijkstra in the battery charging case, the latency was significantly better. The delivery rate was better for the triangles drone network (Figures 13 and 14), with the exception of the Spray and Wait algorithm in the battery changing case. For PRoPHET, Max Prop, and TD-Drone Dijkstra in the battery charging case, the latency was significantly better.

> Latency in the case of changing the battery in the triangles drone network was generally better (Figure 15), with the exception of Spray and Wait, and in the case of Max Delivery a significantly better latency was obtained. In the case of battery charging (Figure 16), latency was better in the triangles drone network for three algorithms (Epidemic, Max Delivery, and TD-Drone Dijkstra), and, for the other three algorithms, the latency was better in the squares drone network. Latency in the case of changing the battery in the triangles drone network was generally better (Figure 15), with the exception of Spray and Wait, and in the case of Max Delivery a significantly better latency was obtained. In the case of battery charging (Figure 16), latency was better in the triangles drone network for three algorithms (Epidemic, Max Delivery, and TD-Drone Dijkstra), and, for the other three algorithms, the latency was better in the squares drone network.

**Figure 13.** Delivery rate, battery changing, squares vs. triangles. **Figure 13.** Delivery rate, battery changing, squares vs. triangles.

Squares Triangles

**Figure 14.** Delivery rate, battery charging, squares vs. triangles.

0

0.2

0.4

Delivery rate

0.6

0.8

**Figure 13.** Delivery rate, battery changing, squares vs. triangles.

liver only on the shorter distances due to buffer restrictions. It is known that Dijkstra's

The delivery rate was better for the triangles drone network (Figures 13 and 14), with the exception of the Spray and Wait algorithm in the battery changing case. For PRoPHET, Max Prop, and TD-Drone Dijkstra in the battery charging case, the latency

Latency in the case of changing the battery in the triangles drone network was generally better (Figure 15), with the exception of Spray and Wait, and in the case of Max Delivery a significantly better latency was obtained. In the case of battery charging (Figure 16), latency was better in the triangles drone network for three algorithms (Epidemic, Max Delivery, and TD-Drone Dijkstra), and, for the other three algorithms, the laten-

algorithm results in the shortest path, and, therefore, the best latency.

was significantly better.

0 0.2 0.4 0.6 0.8 1

Delivery rate

cy was better in the squares drone network.

Squares Triangles

**Figure 14.** Delivery rate, battery charging, squares vs. triangles. **Figure 14.** Delivery rate, battery charging, squares vs. triangles.

**Figure 15.** Average latency: battery changing, squares vs. triangles. **Figure 15.** Average latency: battery changing, squares vs. triangles. **Figure 15.** Average latency: battery changing, squares vs. triangles.

**Figure 16.** Average latency: battery charging, squares vs. triangles. **Figure 16.** Average latency: battery charging, squares vs. triangles.

**Figure 16.** Average latency: battery charging, squares vs. triangles. In essence, the delivery rate was considerably better for the triangles drone networks, and latency was generally better for the triangles case. The significant disad-In essence, the delivery rate was considerably better for the triangles drone networks, and latency was generally better for the triangles case. The significant disad-In essence, the delivery rate was considerably better for the triangles drone networks, and latency was generally better for the triangles case. The significant disadvantage of the triangles drone network is that double the number of drones is needed to cover

vantage of the triangles drone network is that double the number of drones is needed to

vantage of the triangles drone network is that double the number of drones is needed to cover approximately the same surface, and double the number of loading/changing sta-

Epidemic, Spray and Wait, PRoPHET, and MaxProp are classic algorithms used for DTN. However, in our case, as shown, the time-dependent Dijkstra algorithm adaptation can be successfully used because since the flight timetables are known. There are numerous advantages of the TD-Drone Dijkstra algorithm: an exact and optimum route is a priori calculated, ensuring the fastest time of delivery from departure to destination if the route exists; a message is not unnecessarily sent in the network if no route exists from departure to destination; multiple copies of the messages are not unnecessarily spread through drone and station buffers, resulting in unnecessary overloading of the buffers; and, finally, the rate of delivery success is maximized. The drawback of the TD-Drone Dijkstra algorithm is that the route is calculated using the information about the operating drones and stations at the moment of route calculation and, if a drone or a station from the route is down on this route, the message does not reach the destination. All

Epidemic, Spray and Wait, PRoPHET, and MaxProp are classic algorithms used for DTN. However, in our case, as shown, the time-dependent Dijkstra algorithm adaptation can be successfully used because since the flight timetables are known. There are numerous advantages of the TD-Drone Dijkstra algorithm: an exact and optimum route is a priori calculated, ensuring the fastest time of delivery from departure to destination if the route exists; a message is not unnecessarily sent in the network if no route exists from departure to destination; multiple copies of the messages are not unnecessarily spread through drone and station buffers, resulting in unnecessary overloading of the buffers; and, finally, the rate of delivery success is maximized. The drawback of the TD-Drone Dijkstra algorithm is that the route is calculated using the information about the operating drones and stations at the moment of route calculation and, if a drone or a station from the route is down on this route, the message does not reach the destination. All

(63 vs. 65).

(63 vs. 65).

approximately the same surface, and double the number of loading/changing stations is required, although the total number of fixed communication stations is similar (63 vs. 65).

Epidemic, Spray and Wait, PRoPHET, and MaxProp are classic algorithms used for DTN. However, in our case, as shown, the time-dependent Dijkstra algorithm adaptation can be successfully used because since the flight timetables are known. There are numerous advantages of the TD-Drone Dijkstra algorithm: an exact and optimum route is a priori calculated, ensuring the fastest time of delivery from departure to destination if the route exists; a message is not unnecessarily sent in the network if no route exists from departure to destination; multiple copies of the messages are not unnecessarily spread through drone and station buffers, resulting in unnecessary overloading of the buffers; and, finally, the rate of delivery success is maximized. The drawback of the TD-Drone Dijkstra algorithm is that the route is calculated using the information about the operating drones and stations at the moment of route calculation and, if a drone or a station from the route is down on this route, the message does not reach the destination. All the routes passing through the station or drone that is down are compromised until the fault is detected. Moreover, this problem reappears when the drone/station is fixed until the moment this information is updated. However, the chance of this problem occurring is low and, if it appears, it may be fixed in time following the repair of the drone/station or when the current status of the network is updated. Using the example of the Spray and Wait or Epidemic algorithms presented in this paper, any message has the chance to reach the destination even if drones or stations are down because a copy of the message is spread in the network.

#### **4. Conclusions and Future Work**

This paper presents a novel method of communication in quarantined or isolated areas, or areas with technical damage, using networks of drones that fly based on a wellestablished mission plan and schedule, on 2D surfaces covered by squares. Two situations of drone battery management—charging and battery changing stations—were investigated.

A network of square cells with two drones in each cell was proposed to cover a geographical area. The drone network was simulated based on input data from experimental flight tests of a quadcopter using six routing algorithms.

A TD-Drone Dijkstra algorithm (single-copy algorithm) and multiple-copies algorithms were proposed to simulate a Delay Tolerant Network of drones. Results showed a delivery rate ranging from 0.166 to 0.954 in the drone network with battery changing, and from 0.135 to 0.540 in the drone network with battery charging. The best latency of 0.43 h for a drone network with battery changing was obtained using the TD-Drone Dijkstra algorithm. Thus, the best results were obtained for the TD-Drone Dijkstra algorithm, which was able to deliver most of the data packages in the shortest time.

The traditional DTN algorithms, such as Epidemic, Spray and Wait, and MaxProp, produced lower results due to the small number of contacts between nodes, and a low number of message exchanges. The fastest communication was established for the drone squares network with battery changing. It was found that the battery change scenario led to an increase in the delivery rate of 76% compared to the battery charge scenario.

It was found that the battery change scenario led to an increase in the delivery rate of ~200% compared to the battery charge scenario.

The fastest communication was found for the drone triangular network with battery charging, and for the drone square network with battery changing. However, the drone square network is considerably cheaper than the drone triangular network. Thus, if cost is not an issue, the triangle network of drones may be implemented because better performances can be achieved. However, if there is a budget constraint, then the square network type is more suitable.

In future work, we will aim to apply this network of drones to parcel delivery during emergencies in remote quarantined zones.

**Author Contributions:** Conceptualization, A.M.D. and R.U.; methodology, R.U. and A.M.D.; software, A.M.D. and C.-¸S.N.; validation, R.U. and A.M.D.; formal analysis, A.M.D. and R.U.; investigation, R.U. and A.M.D.; resources, R.U., C.-¸S.N. and A.M.D.; data curation, C.-¸S.N. and A.M.D.; writing—original draft preparation, R.U., A.M.D. and C.-¸S.N.; writing—review and editing, R.U. and A.M.D.; visualization, A.M.D. and R.U.; supervision, A.M.D. and R.U.; project administration, R.U. and A.M.D.; funding acquisition, A.M.D., C.-¸S.N. and R.U. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by University Transilvania of Bras, ov.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors acknowledge of Transilvania University of Bra¸sov for providing the infrastructure used in this work.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

