1

Beyond Empirical Models: Pattern Formation
Driven Placement of UAV Base Stations
Jiaxun Lu, Shuo Wan, Xuhong Chen, Zhengchuan Chen, Member, IEEE,

arXiv:1709.03271v1 [] 11 Sep 2017

Pingyi Fan, Senior Member, IEEE, Khaled B. Letaief, Fellow, IEEE

Abstract
This work considers the placement of unmanned aerial vehicle base stations (UAV-BSs) with
criterion of minimum UAV-recall-frequency (UAV-RF), indicating the energy efficiency of mobile UAVs
networks. Several different power consumptions, including signal transmit power, on-board circuit power
and the power for UAVs mobility, and the ground user density are taken into account. Instead of
conventional empirical stochastic models, this paper utilizes a pattern formation system to track the
instable and non-ergodic time-varying nature of user density. We show that for a single time-slot,
the optimal placement is achieved when the transmit power of UAV-BSs equals their on-board circuit
power. Then, for multiple time-slot duration, we prove that the optimal placement updating problem is a
nonlinear dynamic programming coupled with an integer linear programming. Since the original problem
is NP-hard and can not be solved with conventional recursive methods, we propose a sequential-Markovgreedy-decision method to achieve near minimal UAV-RF in polynomial time. Further, we prove that the
increment of UAV-RF caused by inaccurate predicted user density is proportional to the generalization
error of learned patterns. Here, in regions with large area, high-rise buildings or low user density, large
sample sets are required for effective pattern formation.
Index Terms
Aerial base-station, air-to-ground communication, time-varying user density, pattern formation,
Pareto-optimality, sampling number, UAV deployment.
Jiaxun Lu, Shuo Wan, Xuhong Chen and Pingyi Fan are with Tsinghua National Laboratory for Information Science and
Technology(TNList) and the Department of Electronic Engineering, Tsinghua University, Beijing, P.R. China, 100084. E-mail:
{lujx14, s-wan13, chenxh13}@mails.tsinghua.edu.cn, fpy@tsinghua.edu.cn.
Zhengchuan Chen is with the College of Communication Engineering, Chongqing University, Chongqing, P.R. China, 400044.
E-mail: chenzc.ee@gmail.com.
Khaled B. Letaief is with the Department of Electrical and Computer Engineering, HKUST, Hong Kong (email:
eekhaled@ust.hk) and Hamad bin Khalifa University, Qatar (email: kletaief@hbku.edu.qa).

2

I. I NTRODUCTION
Unmanned aerial vehicle base stations (UAV-BSs) have been considered as a promising
solution to provide wireless coverage in a rapid manner. In this system, UAVs are often powered
by batteries [1], which limits its life-time. Here, UAV-BSs are usually recalled periodically and
there is an urgent need to increase the energy-efficiency of UAVs systems. A general strategy to
improve energy-efficiency is adjusting the placement of UAVs according to ground user density,
working environment and desired transmit data rate etc. [2]–[4].
In the literatures, the placement of UAV-BSs has been studied under a broad range of aspects
[2]–[8]. The optimal hovering altitude that maximizes coverage radius was discussed in [5].
Later, the authors in [2], [8] investigated the optimal hovering altitude and coverage radius.
In [6], the authors analyzed the optimal placement of UAV-BSs simultaneously maximizing the
number of covered users and energy-efficiency. In [3], [4] and [7], they discussed the relationships
among the optimal placement, minimum required number of UAV-BSs, as well as the density
of ground users. The main focus of these works is on minimizing transmit power. In fact, the
on-board circuit power consumption, related to rotors, computational chips and gyroscopes etc.,
and the potential mobility power consumption of UAVs may also affect the life-time of network
and should be taken into account in the viewpoint of system. Up to now, the works on energyefficient placement of UAV-BSs considering on-board circuit power and mobility power are quite
limited. A power model investigating the peculiar features of UAVs, like available energy, weight,
maximum speed, etc. was formulated in [9]. The on-board circuit power and mobility power of
UAVs were addressed in [10] and [11], separately. These works provide comprehensive analyses
based on a time-invariant density of ground users. However, in practice, due to the directional
shift of people consequences for transportation design or epidemic control etc., the human flow
usually follows Turing pattern [12], which leads to time-varying, instable and non-ergodic density
of ground users, raising some new challenges to the placement of UAV-BSs.
In conventional stable and ergodic cases, one can develop an empirical stochastic model
for density of ground users (e.g., Markov modulated) where the statistical parameters can be
estimated from data [13]. Such models may not be applied to the instable scenarios. Fortunately,
the developing data-driven methods cast a bright light on challenges in wireless communication.
For example, useful information extraction with Data-mining approach was utilized in [14] to
enhance the caching performance of wireless networks. Machine-learning techniques have been

3

applied to decision-making and feature classification in cognitive radio problems (see [15] and
the references therein). Some tutorial works on applying data-driven methods to the wireless
domain were presented in [16] and [17]. These works deliver some useful motivations to us.
Recently, it was shown in [18] and [19] that using machine learning, it is possible to construct a
more intelligent context-aware pattern formation system by predicting future situations as well
as monitoring the current state.

Statistics

Parameters Empirical
Model

Data
Pattern
Formation

Pattern

Decision
Model

Output
Placement
sets

Fig. 1. Pattern Formation aided framework of our considered system.

As shown in Fig. 1, we propose a framework, using a pattern formation module rather than a
statistical empirical model, to track the instable and non-ergodic time-varying nature of ground
user densities. Then, we address the importance of accurate pattern formation and consider the
decision model on optimal placement that maximizes the life-time of the mobile UAVs network.
In this case, the on-board circuit power and the potential mobility power of UAVs are also
considered. This paper focuses on the downlink of UAV-BSs, in which each of the ground users
is served with fixed data rate. In the system, the considered duration is partitioned into continuous
time-slots. At the beginning of each time-slot, UAV-BSs are allowed to decide whether or not to
update their placement according to ground user density. To our best knowledge, this paper is
one of the first comprehensive studies on the joint optimal deployment of UAV-BSs and pattern
formation in scenario with instable and non-ergodic time-varying density of ground users.
To characterize the life-time of mobile UAVs network, we employ the notion of UAV-recallfrequency (UAV-RF), the frequency of the active UAVs run out of batteries, as the physical
index. That is, maximizing life-time is equivalent to minimizing UAV-RF. In this direction, we
first consider the optimal placement of UAV-BSs that minimizes UAV-RF in one time-slot and
then extend the discussion to multi-slot duration. In fact, the UAV-RF in one time-slot can
be treated as static UAV-RF. In this case, by analyzing the coverage scenario with one single
UAV, we prove that the optimal hovering altitude minimizing transmit power is proportional to

4

the coverage radius, and the slope is only determined by communication environment (highrise urban, dense urban and urban, etc.), which is a general extension of previous results in
[5] by considering the density of users inside the coverage of UAV-BSs. More specifically, in
environment with high-rise buildings, the slope is large, and hence UAVs are supposed to fly
higher compared with environment with low-rise buildings.
By applying the derived optimal hovering altitude, we investigate the static UAV-RF versus environment, coverage parameters and on-board circuit power, where coverage parameters include
the coverage radius, user density and desired data rate. Analytical results demonstrate that: 1)
The minimal static UAV-RF is achieved when transmit power equals the on-board circuit power;
2) The minimal static UAV-RF becomes large in scenarios with high-rise buildings, high onboard circuit power, and large user density and data rate. This indicates that limiting on-board
circuit power can effectively prolong the life-time of mobile UAVs network. Compared with the
optimal coverage radius given in [8], our results provide more insights on the design of UAVs
networks by investigating the on-board circuit power.
For the multiple time-slot case, it requires to decide the optimal placement updating epochs of
UAV-BSs in cases with instable time-varying density of ground users. The corresponding UAVRF is denoted as dynamic UAV-RF, which is relevant to transmit power, on-board circuit power
and potential mobility power of UAV-BSs. We show that the placement optimization problem is
a multi-stage decision process and can be written as a nonlinear dynamic programming (NLDP)
coupled with an inherent integer linear programming (ILP). For the inherent ILP, we propose
a polynomial time solution by transforming it into a standard assignment problem with some
amendments. However, for the NLDP, we show that it can not be solved with conventional
methods, such as recursive manners, because the number of update epochs is unknown in
advance.
Noticing the NP-hardness of NLDP, we shall design a sequential-Markov-greedy-decision (SMGD) method to find the near-optimal solution by utilizing the notion of Pareto-optimality,
which is proved to be with polynomial complexity. Simulation results show that our proposed
S-MGD method can stably achieve near-optimal performance in terms of minimum dynamic
UAV-RF. In particular, when mobility power is extremely low compared with on-board circuit
power and transmit power, the proposed S-MGD method updates the placement of UAV-BSs
at the beginning of each time-slots. By contrast, in cases with extremely high mobility power,
UAV-BSs hold their placement during the considered duration.

5

Finally, the relationships among sampling number, density pattern accuracy and increment of
UAV-RF are characterized in detail. These results imply that in subregions with large area, highrise buildings and low user density, large sample sets are needed for effective pattern formation
and reducing UAV-RF. Specifically, we first prove that the increment of UAV-RF caused by
inaccurate density patterns is proportional to the generalization error. Then, we theoretically
derive the minimum sampling number of each subregion with the Vapnik-Chervonenkis (VC)
theorem [20], where the overall increased UAV-RF is upper-bounded.
The rest of this paper is organized as follows. In Section II, the system model is introduced,
and an available density pattern is presented as the test set in this paper. Then, the optimal
placement of UAV-BSs minimizing static UAV-RF is discussed in Section III. In Section IV, the
S-MGD based placement updating method is presented to minimize dynamic UAV-RF. Section
V analyzes the effects of the accuracy of pattern formation system. In Section VI, the validity
of previous theoretical results and the effectiveness of our proposed methods are verified by
numerical results. Finally, conclusions are given in Section VII.
II. C OVERAGE M ODEL , A IR - TO -G ROUND C HANNEL

AND

U SER D ENSITY PATTERN

In this section, we shall first illustrate the downlink coverage model of UAV-BSs, where
instable time-varying user density is considered. Then, we introduce the air-to-ground (A2G)
channel and show the existence of optimal hovering altitude of UAV-BSs. Finally, an available
pattern of previous time-varying user density is presented as the test set in this paper.
A. UAV-BS Coverage Model
As shown in Fig. 2, a geographical area is divided into several subregions according to
the municipal planning of modern cities and their different ground user density patterns [21].
Consider a length T duration, let the time-varying density of ground users be λβ (t), where
β ∈ {1, 2, · · · , κ} and t ∈ [0, T ] index subregion and time, respectively. κ is the number of
subregions. For the notation simplicity, unless otherwise stated, we drop the time index from the
equations. To balance the number of active UAV-BSs between adjacent time-slots, a recall and
supplement center (RSC) is deployed in considered area.
We focus on the downlink of network in which UAV-BSs adopt a frequency division multiple
access technique to serve each of the ground users with fixed data rate C. UAV-BSs assign
individual frequency bands to mobile ground users, and hence the frequency interference between

6

UAV

hb

E
R

q

r

C

Rb
O

T
C
R

Comprehensive

Resident

T

E

Entertainment

Transport

O

Office

Fig. 2. A typical area classified into different subregions according to user density patterns and the UAV coverage model therein.

UAV-BSs is avoided. We assume that the transmit power of each UAV-BS and the available
bandwidth are sufficient to meet the rate requirement of users. In this case, we consider a diskcovering model where UAV-BSs in the same subregion provide equal coverage radius [7], [22].
Here, considering the fact that the area of overlaps between adjacent disks are proportional to
the area of disks, without loss of generality, we express the required number of UAV-BSs in the
β-th subregion at t as1
Nβ (t) =

Sβ
,
πRβ2 (t)

(1)

where Sβ is the area of the β-th subregion and Rβ (t) is the corresponding coverage radius.
B. Air-to-Ground Channel
The A2G channel can be characterized into line-of-sight (LOS) link or non-line-of-sight
(NLOS) link [5], and the path loss therein can be given by



(4πf /c)2 r 2 + h2 η0 ,
β
Lξ (r, hβ ) =


(4πf /c)2 r 2 + h2 η1 ,
β

ξ=0

(2)

ξ = 1,

where ξ = 0 and ξ = 1 denote LOS link and NLOS link, respectively. f is the carrier frequency
and c is the traveling speed of light. r ∈ [0, Rβ ] is the distance between the user of interest and
the projection of UAV-BS on ground, and hβ is the hovering altitude of UAV-BSs. η0 and η1 are
the excessive path loss on the top of the free space path loss (FSPL) for LOS and NLOS links,
1

Minimizing required number of UAV-BSs can be formulated as the geometric disk cover problem [23], and can be optimally

solved by core-sets method [24].

7

determined by communication environment (suburban, urban, dense urban, high-rise urban or
others). Typically, η1 ≫ η0 as the obstacles in propagation paths greatly enhance the path loss
of NLOS link.
The average path loss of A2G channel is also determined by the LOS probability [5]
P0 (r, hβ ) =

1
,
1 + aexp(−b[θ − a])

where a and b are constants determined by environment, and θ =

180
π

(3)
tan−1 (hβ /r) is the elevation

angle shown in Fig. 2. Then, 1 − P0 (r, h) denotes the NLOS probability. The average path loss
of A2G channel can be derived as
L̄(r, hβ ) = P0 (r, hβ )L0 (r, hβ ) + (1 − P0 (r, hβ )) L1 (r, hβ )


= (4πf /c)2 r 2 + h2β η1 + P0 (r, hβ )(η0 − η1 ) .
{z
}
{z
}|
|
FSPL

(4)

average excessive path loss

This clearly characterizes the individual effects of FSPL and average excessive path loss. The
first part accounts for FSPL, which monotonically increases with hβ due to the growing distance
between UAV and user; the second part, which represents the average excessive path loss, is
monotonically decreasing with hβ . This is because large hβ leads to high LOS probability of
A2G channel. For a specific coverage radius Rβ , (4) implies the existence of optimal hovering
altitude minimizing average path loss of A2G channel. Besides, to minimize the UAV-RF of
considered area, the optimal coverage radius also needs to be jointly considered.
C. An Available Density Pattern as the Test Set
The empirical average traffic amount of the β-th subregion can be reconstructed as {xrβ [n]}
by the inverse discrete Fourier transformation (IDFT) [21]:
xrβ [n] =

N −1
γr,β X r
X [k] exp(2πjkn/N),
N k=0 β

(5)

where n ∈ [0, N] is the sampling index, and with the sampling period µ = 10 min, the sampling
number in 4 weeks is N = 4032. γr,β is the reconstruction scaling factor at the β-th subregion,
and Xβr [k] is the IDFT coefficient expressed by


Xβ [k], k = 0, 4, 28, 56, N − 4, N − 28, N − 56
r
Xβ [k] =

0,
otherwise.

(6)

8

Xβ [k] is the frequency-domain coefficients of DFT(xβ [n]), where xβ [n] is the sampled timedomain traffic amount in the β-th subregion. According to the properties of IDFT, Xβ [k] =
Xβ† [N − k] for the real time-domain traffic amount, where (·)† denotes the conjugate transpose
of (·). Hence, only parts of the coefficients and the scaling factor γr,β are listed in Table I.
TABLE I
C OEFFICIENTS OF R ECONSTRUCTION F UNCTION

Subregions

γr (×1011 )

k = 0 (×10−4 )

k=4

k = 28

k = 56

E

8.35

3.24

0.06e−0.3j

0.5e2.36j

0.08e0.69j

R

17.4

2.73

0.04e−1.02j

0.27e1.72j

0.17e1.35j

T

4.32

3.73

0.1e1.04j

0.38e2.53j

0.28e2.46j

O

5.23

4.63

0.21e1.21j

0.56e2.52j

0.2e0.29j

C

17.4

2.85

0.04e0.35j

0.3e2.19j

0.15e1.11j

Then, with the assumption that each of the ground users is served with fixed data rate C, the
average user density in the β-th subregion at t can be immediately expressed as
λβ (t) =

1
xrβ [⌊t/µ⌋] .
CSbs

(7)

Sbs is the area of the coverage of base stations shown in [21]. The normalized density patterns in
one week corresponding to five considered subregions are depicted in Fig. 3. It can be observed
that the density patterns are instable and show the breath-out mode in the considered duration. In
this case, to provide optimal coverage in each subregion, the placement of UAV-BSs is supposed
to be updated with respect to the learned time-varying ground user densities. However, updating
the placement requires the mobility of UAV-BSs, which may cost much energy and reduce the
life-time of mobile UAVs network [9]. Therefore, to minimize the UAV-RF of network, the
optimal placement of UAV-BSs in one time-slot and the placement updating strategy in the
considered duration should be jointly optimized.
III. O PTIMAL P LACEMENT

OF

UAV-BS S

IN SINGLE TIME - SLOT

In this section, we first formulate the optimal UAV-BSs placement problem in single time-slot.
Then, we show that the original problem can be solved by separately investigating the optimal
hovering altitude and optimal coverage radius. Finally, we present the optimal condition with
respect to on-board circuit power.

9

14:00-19:00 Monday

E

R

T

O

C

0

1

2

3

4

5

6

7

Time / day

Fig. 3. Density patterns of ground users in different subregions.

A. Problem Formulation
Recall that the downlink rate is fixed to C. Denote the allocated transmit power to the interested

user located at r be Ptr,user,ξ r, hβ . Then, following Shannon formula, one has
!
Ptr,user,ξ r, hβ

,
(8)
C = W log2 1 +
Lξ r, hβ N0 W
where N0 is the noise power spectrum density, and W is the allocated bandwidth to interested

user. The relevant transmit power can be immediately derived as


Ptr,user,ξ r, hβ = Lξ (r, hβ )N0 W 2C/W − 1 .

(9)



P̄tr,user r, hβ ) = L̄(r, hβ N0 W 2C/W − 1 .

(10)

Accordingly, the average transmit power is given by

The expected value of transmit power of UAV-BS is the integral of average transmit power
relevant to all users inside the coverage. That is,
Z Rβ


2πr P̄tr,user r, hβ dr.
Ptr Rβ , λβ , hβ = λβ

(11)

0

Let the on-board circuit power and the battery capacity of one single UAV-BS be Pcu and Eb ,

respectively. Then, the average life-time of UAV-BSs in the β-th subregion is

Eb / Ptr (Rβ , λβ , hβ ) + Pcu ,

(12)

10

and accordingly the static UAV-RF at the β-th subregion can be expressed as
Φst,β (t) =

Nβ (t) (Ptr (Rβ , λβ , hβ ) + Pcu )
.
Eb

(13)

From (13), it can be observed that static UAV-RF is determined by the number of UAV-BSs
and the power relevant to signal transmission and on-board circuit. Besides, the total consumed
power of the mobile UAVs network can be written as


κ
κ
X
X

Φst,β (t).
Nβ (t) Ptr Rβ , λβ , hβ + Pcu = Eb

(14)

β=1

β=1

That is, the total UAV-RF of our considered area takes both the power consumed by one single
UAV-BS and the number of UAV-BSs into account. Compared with the total consumed power,
the notion of UAV-RF characterizes the frequency of active UAVs run out of batteries and is a
more comprehensive indicator of the life-time of mobile UAVs network.
As illustrated in Section II-B, to minimize static UAV-RF of considered area in single time-slot,

the altitudes and coverage radii of UAV-BSs need to be jointly considered. That is,
P1 :

min
{Rβ , hβ }

κ
X

Φst,β (t)

(15)

β=1

s.t. hβ ≥ 0, Rβ > 0.

(16)

In P1 , (16) corresponds to the non-negativity of hovering altitudes and coverage radii of UAVBSs. According to (13), it can be observed that the optimization variables {Rβ , hβ } are coupled,
and hence directly solving problem P1 is difficult. In the following, we consider a two-stage
method instead.
B. Optimal Hovering Altitude of UAV-BSs
We firstly consider the optimal hovering altitude minimizing the transmit power with a fixed
coverage radius, which can be expressed as
P1-A : min Ptr Rβ , λβ , hβ
{h β }
s.t. hβ ≥ 0.



(17)
(18)

(18) reflects the non-negativity of hovering altitude. Then, the solution to P1-A can be summarized
as follows.

11

Lemma 1. 1) The transmit power of single UAV-BS can be expressed as


Ptr Rβ , λβ , hβ = γtr Ptr,1 hβ /Rβ ,

Ptr,1 (hβ /Rβ ) = Ptr 1, 1, hβ /Rβ

(19)



is the transmit power corresponding to hovering altitude

hβ /Rβ , when Rβ , λβ and C are normalized. γtr = λβ Rβ4 2C/W − 1 is the relevant scaling
factor.

2) The optimal hovering altitude corresponding to specific Rβ is h∗ = Rβ h∗β,1 , where h∗β,1 is
the optimal hovering altitude that minimizes Ptr,1 (hβ,1 ) and is only determined by communication
environment.
Proof. With (11), we can derive that
Z

Ptr,1 (hβ /Rβ ) =

1

2πr L̄(r, hβ /Rβ )N0 W dr.

(20)

0

Besides, (4) can be rewritten as L̄(r, hβ ) = Rβ2 L̄(r/Rβ , hβ /Rβ ). The transmit power of one
single UAV-BS can be expressed as
Ptr (Rβ , λβ , hβ ) = λβ

Z

Rβ

2πr P̄u(r, hβ ) dr

0

=

λβ Rβ2

2

C/W

−1



Z

Rβ

2πr L̄(r/Rβ , hβ /Rβ )dr

0

= λβ (β, t)Rβ4 2C/W − 1



Z

(21)

1

2πr L̄(r, hβ /Rβ )dr

0


= λβ Rβ4 2C/W − 1 Ptr (1, 1, hβ /Rβ ).

This completes the proof of (19). Then, one can immediately have

∂Ptr Rβ , λβ , hβ
∂Ptr,1 (hβ /Rβ )
=0⇔
= 0.
∂hβ
∂hβ

(22)

Noticing that Ptr,1 (hβ /Rβ ) only accounts for the communication environment, h∗β = Rβ h∗β,1
naturally follows.
Lemma 1 clearly indicates the individual effects of communication environment and coverage
parameters on optimal hovering altitude. That is, h∗β,1 accounts for the environmental statistics,
while γtr,1 explains the effects of coverage parameters. This implies that solving P1-A is equivalent
to finding the optimal hovering altitude that minimizes Ptr,1 (hβ,1 ). That is,


∂Ptr,1 (hβ,1 )
∗
hβ,1 = arg
=0 .
∂hβ,1
hβ,1

(23)

12

From (20), we have that
∂Ptr,1 (hβ,1 )
=0⇔
∂hβ,1

Z 1n
0


2hβ,1 η1 + P0 (r, hβ,1 )(η0 − η1 ) +
2

r +

h2β,1

where based on (3),



o
∂P0 (r, hβ,1 )
(η0 − η1 ) rdr = 0,
η1 +
∂hβ,1

180brP0 (r, hβ,1 )
∂P0 (r, hβ,1 )
=
(1 − P0 (r, hβ,1 )).
∂hβ,1
π(r 2 + h2β,1 )

(24)

(25)

Substituting (24) and (25) into (23), one can get h∗β,1 . However, it is overwhelming to obtain
explicit solution of (23). As an alternative, we propose a binary search algorithm to calculate
the optimal hovering altitude, as shown in Algorithm 1. Note that the precision ǫ = 10−3 and
the iteration scaling factor γi = 10 shown at line 0 of Algorithm 1 can be readily replaced with
any other values that satisfy ǫ > 0 and γi > 1.
Algorithm 1 Optimal Hovering Altitude
0: Initialize
Environmental parameters: η1 , η2 , a and b;
Input coverage parameters: Rβ , λβ , C;
Initialize iteration parameters: h∗β,1 = hmin = 0, hmax = 1;

2:

Set the precision ǫ = 10−3 and iteration scaling factor γi = 10.


∂Ptr,1 (hβ )
∂Ptr,1 (hβ )
|
|
while
hβ =hmax
hβ =hmin ≥ 0 do
∂hβ
∂hβ

3:

end while

4:

while

1:

hmax = γi hmax ;

∂Ptr,1 (hβ )
|hβ =h∗β,1
∂hβ

≥ ǫ do

5:

h∗β,1 = (hmax + hmin ) /2;

6:

if

7:
8:
9:

∂Ptr,1 (hβ )
|hβ =h∗β,1
∂hβ

≥ 0 then

hmax = h∗β,1 ;
else
hmin = h∗β,1 ;

10:

end if

11:

end while
Output

The optimal hovering altitude h∗β = Rβ h∗β,1 .

13

C. Optimal Coverage Radius of UAV-BSs
According to Lemma 1, the optimal hovering altitudes of UAV-BSs are determined by the
coverage radius. Then, P1 can be immediately rewritten as
κ
X
Φst,β (t)
P1-B : min
{Rβ } β=1

(26)

s.t. h∗β = Rβ h∗β,1 ,

(26a)

Rβ > 0.

(26b)

(26a) is the optimal hovering altitude corresponding to Rβ , and (26b) reflects the non-negativity
of coverage radius. As shown in (26), the static UAV-RF of considered area is the sum of
individual static UAV-RF in each subregion. Hence, solving P1-B is equivalent to minimizing
Φst,β (t) for all β ∈ {1, 2, · · · , κ}, separately.
Substituting (1) and (19) into (13), the static UAV-RF can be expressed as



Sβ Pcu
C/W
∗
2
Φst,β (t) =
+ λβ 2
− 1 Ptr,1 (hβ,1 )Rβ .
πEb Rβ2

(27)

Thus, with the theoretical results shown in Section III-B, the optimal placement of UAV-BSs
that minimizes static UAV-RF in the considered area can be summarized as follows.
Theorem 1. 1) The optimal coverage radius and hovering altitude of UAV-BS in the β-th
subregion at t is
Rβ∗ (t)

=

s
4

λβ

(t) (2C/W

Pcu
− 1) Ptr,1 (h∗β,1 (t))

(28)

and
h∗β (t) = Rβ∗ (t)h∗β,1 (t),

(29)

respectively. h∗β,1 (t) is given by (23).
2) The optimal placement of UAV-BSs in the β-th subregion is achieved when the on-board
circuit power of one UAV-BS equals its transmit power. That is,

where γtr∗ (t) = λβ (t)Rβ∗4 (t) 2C/W

Pcu = γtr∗ (t)Ptr (h∗β,1 (t)),

−1 .

(30)

Proof. Denote the minimal static UAV-RF as Φ∗st,β (t). Since both the on-board circuit power and
user density in (27) are positive, one has
2Sβ q
Φ∗st,β (t) =
λβ (2C/W − 1) Pcu Ptr,1 (h∗β,1 ),
πEb

(31)

14

and Φ∗st,β (t) is achieved when Rβ (t) = Rβ∗ (t), which is shown in (28). According to Lemma 1,
the optimal hovering altitude corresponding to Rβ∗ (t) is given by (29). Substitute (28) into (19),
(30) can be easily proved.
Theorem 1 not only provides the optimal hovering altitude and the optimal coverage radius,
but also points out the optimal transmit power and circuit power allocation. The former result
provides valuable insights for UAV-BSs positioning while the latter result presents an efficient
resource allocation for multi-function usage of UAV-BSs. The physical meaning of Theorem 1 is
intelligible. When on-board circuit power is high, large coverage radius can decrease the number
of active UAV-BSs. According to (13), small Nβ (t) decreases the effects of high Pcu on static
UAV-RF, and hence Φst,β (t) is reduced. By contrast, when on-board circuit power is low, small
coverage radius can decrease transmit power, which also decreases Φst,β (t). Specifically, when
Pcu = 0, we have Rβ∗ = 0 and h∗β = 0. That is, users can connect to UAV-BSs just at their
balance positions. In this way, enlarging the number of UAVs does not consume any on-board
circuit power, while the transmit power is saved.
IV. DYNAMIC P LACEMENT

OF

UAV-BS S

IN

C ONSIDERED D URATION

This section shall consider the placement of UAV-BSs in time-dimension. As illustrated in
Section II-C, there exists a trade-off between updating the placement of UAV-BSs and reducing
the dynamic UAV-RF. In this case, we firstly formulate the optimal placement updating problem
in terms of minimum dynamic UAV-RF. We demonstrate that the original problem is NP-hard
and can not be solved in conventional manners. Finally, a sequential method is proposed to
update the UAV-BS placement near-optimally in polynomial time.
A. Formulation of the Optimal Placement Updating
Recall that only at beginning of time-slots, the placement of UAV-BSs can be updated to be
optimal. Denote the update epochs of UAV-BSs in the considered area as τi (i = 0, · · · , Nτ ). Nτ
is the number of update epochs and 0 ≤ Nτ ≤ T /µ. Specifically, let τ0 = 0. Denote the updated
coverage radius of UAV-BSs in the β-th subregion at τi as Rβ (τi ). The corresponding number
of UAV-BSs can be expressed as Nβ (τi ) = Sβ /Rβ (τi ). Furthermore, if Nτ ≥ 1, the number of
UAV-BSs need re-positioning at τi can be given by
ζτ |ττii−1 = max (N(τi−1 ), N(τi )) ,

(32)

15

where N(τi ) =

Pκ

β=1

Nβ (τi ) and i = 1, · · · , Nτ .

Then, the mobility energy of UAV-BSs at τi (i = 1, · · · , Nτ ) can be expressed as [9]
τ

ζτ |τii−1

Ωm (τi ) =

X
l=1

Ph

 ∆hβ (τi , l)
dβ (τi , l)
∆hβ (τi , l)
+I(∆hβ (τi , l))Pa
− 1−I(∆hβ (τi , l)) Pd
, (33)
vh
va
vd

where Ph and vh are the mobility power and velocity for one single UAV-BS in the horizontal
direction, respectively. dβ (τi , l) and ∆hβ (τi , l) are the horizontal moving distance and the variation of the height of the l-th UAV-BS in the β-th subregion at τi , respectively. Pa and va denote
the ascending power and ascending velocity of one single UAV-BS, respectively. Similarly, Pd
and vd denote the descending power and descending velocity. I(∆hβ (τi , l)) is the indicative
function, i.e.
I(∆hβ (τi , l)) =



1, ∆hβ (τi , l) ≥ 0

(34)


0, ∆hβ (τi , l) < 0.

Specifically, Ωm (τ0 ) = 0. Note that Ωm (τi ) (i = 1, · · · , Nτ ) is determined by the total moving
distance of UAV-BSs at τi , related to the placement of UAV-BSs at τi−1 .
We consider the optimal updating strategy as {τ0 , · · · , τNτ }. With (33), the average dynamic
UAV-RF in considered duration can be given by
)
(Z
Z TX
Nτ
κ
κ
T X
X
1
1
1
Ωm (τi ) ,
Φst,β (t)dt +
Φdn,β (t)dt =
Φ̄dn |T0 =
T 0 β=1
T
Eb i=0
0 β=1

(35)

where Φdn,β (t) is the instantaneous dynamic UAV-RF at t. The corresponding optimal placement
updating method can be immediately expressed as
P2 :

min

{τ0 ,··· ,τNτ }

Φ̄dn |T0

s.t. 0 ≤ Nτ ≤ T /µ,

(36)
(36a)

τ0 = 0, τNτ < T,

(36b)

if Nτ ≥ 1, τi < τi+1 (i = 0, · · · , Nτ − 1).

(36c)

(36a) shows the bounds of update epochs. (36b) and (36c) denotes the range and the order of
update epochs, respectively.
Problem P2 is a typical multi-step decision process. Because the objective function given
by (33) is nonlinear, P2 can be solved by nonlinear dynamic programming (NLDP) method.
However, because the number of updating times is unknown in advance, the number of recursion

16

ΨR =



(ψR (1), · · · , ψR (|N(τi ) − N(τi+1 )|)) ,



∅,

if N(τi ) 6= N(τi+1 )

(38)

if N(τi ) = N(τi+1 ).

formula can not be determined as well. Hence, P2 is NP-hard and difficult to be directly solved
as the size of problem grows. Hereinafter, to efficiently find the near optimal solution of P2 , we
propose a new sequential method with polynomial computational complexity.
B. Trajectory Planning
Given the update epochs τi (i = 1, · · · , Nτ ), the mobility energy shown in (33) is only
determined by the updating trajectory of the placement of UAV-BSs. Therefore, we firstly analyze
the optimal trajectory planning method that minimizes the mobility energy of UAV-BSs.
Let the 3D position of the k-th UAV-BS at τi−1 be ψ(k, τi−1 ), where k = 1, · · · , N(τi−1 ).
Then, the corresponding position set of UAV-BSs can be expressed as
Ψ(τi−1 ) = (ψ(1, τi−1 ), · · · , ψ(N(τi−1 ), τi−1 )) , i = 1, · · · , Nτ .

(37)

To balance the number of UAV-BSs before and after updating the placement of UAV-BSs at
τi , the redundant UAV-BSs should be recalled to RSC when N(τi−1 ) > N(τi ), whereas the
additional UAV-BSs should be supplemented by RSC when N(τi−1 ) < N(τi ). Let the position
of RSC be ψR (k). The position set of recalled or supplemented UAV-BSs can be expressed as
(38), which is shown at the top of this page. Then, the position sets of UAV-BSs at τi−1 and τi
can be immediately given by

Ψ̂(τi−1 ) =
and
Ψ̂(τi ) =



Ψ(τi−1 ),


(Ψ(τi−1 ), ΨR ) ,



(Ψ(τi ), ΨR ) ,

Ψ(τi ),

if N(τi−1 ) ≥ N(τi )

(39)

if N(τi−1 ) < N(τi )
if N(τi−1 ) ≥ N(τi )

(40)

if N(τi−1 ) < N(τi ),

i = 1, · · · , Nτ . Obviously, the volume of both the position sets Ψ̂(τi−1 ) and Ψ̂(τi ) is ζτ |ττii−1 . Our
goal is to find the optimal mapping between Ψ̂(τi−1 ) and Ψ̂(τi ) that minimizes Ωm (τi ).

17

According to (33), one can see that finding the optimal mapping is a standard integer linear
programming (ILP). In general, this problem can be solved by using standard ILP solution
methods. However, these solutions may not be efficient as the size of the problem grows. Due
to the potential high number of UAV-BSs, a more efficient technique is needed. To this end, we
transform the trajectory planning problem into a standard assignment problem, which can be


τi+1 3
solved with the Hungarian method in polynomial time O (ζτ |τi ) [25].
Define the mobility energy matrix as C, which is a ζτ |ττii−1 ×ζτ |ττii−1 square matrix. Let ψ̂x,y (k, τi )

and ψ̂z (k, τi ) be the horizontal and vertical coordinate of the k-th element of ψ̂(k, τi ), respectively.
Then, the element at the k-th row and the l-th column of C can be immediately written as
C(k, l) =Ph

||ψ̂x,y (l, τi ) − ψ̂x,y (k, τi−1 )||2
+
vh
ψ̂z (l, τi ) − ψ̂z (k, τi−1 )
−
va
 ψ̂z (l, τi ) − ψ̂z (k, τi−1 )
.
1 − I(ψ̂z (l, τi ) − ψ̂z (k, τi−1 )) Pd
vd

I(ψ̂z (l, τi ) − ψ̂z (k, τi−1 ))Pa

(41)

C(k, l) represents the mobility energy when a UAV-BS moves from ψ̂(k, τi−1 ) to ψ̂(l, τi ). Denote
the corresponding ζτ |ττii−1 × ζτ |ττii−1 assignment matrix as Z, where the element Z(k, l) is 1 if the
UAV-BS at ψ̂(k, τi−1 ) is assigned to ψ̂(l, τi ), or 0 otherwise. In this way, the standard assignment
problem on finding the optimal trajectory can be formulated as
τ

τ

X

X

ζτ |τii−1 ζτ |τii−1

P2-A : min
Z

k=1

τ
ζτ |τii−1

s.t.

X

C(k, l)Z(k, l)

(42)

l=1

τ

ζτ |τii−1

Z(k, l) = 1,

k=1

X

Z(k, l) = 1.

(43)

l=1

Constraint (43) guarantees that the assignment between Ψ̂(τi−1 ) and Ψ̂(τi ) is one-to-one.
C. Placement Updating Strategy
Following the optimal trajectory planning method, this part gives the optimal updating strategy
{τ0 , · · · , τNτ } by solving problem P2 . To this end, we first illustrate the Markov nature of P2 ,
as shown in the following Lemma.
Lemma 2. The static UAV-RF at t can be expressed as
Φst,β|(τ0 ,··· ,τi ) (t) = Φst,β|τi (t),

(44)

18

where τi = max{τi |τi ≤ t}, i = 0, · · · , Nτ . Φst,β|(τ0 ,··· ,τi ) (t) and Φst,β|τi (t) are the corresponding
static UAV-RFs at time t when (τ0 , · · · , τi ) and τi are given, respectively.
Proof. Since the placement UAV-BSs has been updated at τi , the coverage radius at time t can
be given by Rβ (t) = Rβ (τi ). In addition, because as shown in (27), the static UAV-RF is only
determined by the coverage radius Rβ (t), (44) follows.
With Lemma 2, we can see that the static UAV-RF at time t is determined by the nearest
update epoch in a Markov mode. This motivates us to investigate the average dynamic UAV-RF
between consecutive update epochs, denoted by
R τi Pκ
β=1 Φst,β (t)dt + Ωm (τi )
τ
Φ̄dn |ττii−1 = i−1
, i = 1, · · · , Nτ .
τi − τi−1

(45)

Considering the fact that the number of update epochs is unknown, we can sequentially decide the
update epochs by minimizing Φ̄dn |ττii−1 . The following Theorem established the Pareto-optimality
of this sequential decision method.
Theorem 2. Given τ0 = 0, the sequential-Markov-greedy-decision (S-MGD) method can be
expressed as
τi = argmin Φ̄dn |ττii−1 , i = 1, · · · , Nτ ,

(46)

τi ∈(τi−1 ,T )

Average dynamic UAV-RF

and the corresponding {τ0 , · · · , τNτ } is a Pareto-optimal placement updating strategy.

12

10

8

6

4
i-1

i

i

+
i

i+1

Time [h]

Fig. 4. A general case of proposed S-MGD method. The optimal update epochs are marked with stars.

Proof. For the convenience of proof, we utilize a general case shown in Fig. 4 to demonstrate
the Pareto-optimality of proposed S-MGD method. The solid line and dotted line illustrate how
Φ̄dn |ττii−1 varies with respect to τi , given τi−1 . According to (46), the optimal update epoch τi is

19

greedily determined and is marked with star in Fig. 4. That is, Φ̄dn |ττii−1 achieves its minimum
τ

at the optimal τi . Similarly, the red dashed line depicts how Φ̄dn |τii+1 varies with τi+1 when τi
τ

is given by (46). Also, the optimal τi+1 that minimizes Φ̄dn |τii+1 is marked with star in Fig. 4.
Then, we prove the Pareto-optimality of Theorem 2. That is, for the average dynamic UAV-RF
between each consecutive update epochs during [0, T ], we can not decrease one of them without
increasing the others.
Take one of the the duration [τi−1 , τi+1 ] (i = 1, · · · , Nτ −1) for example. According to Lemma
τ

2, adjusting τi only affects Φ̄dn |ττii−1 and Φ̄dn |τii+1 . Let τi− and τi+ denote the time belonging to
(τi−1 , τi ) and (τi , τi+1 ), respectively. As shown in Fig. 4, due to τi minimizes Φ̄dn |ττii−1 , one has

−

Φ̄dn |ττii−1 > Φ̄dn |τi
τi−1
(47)
+

τ
Φ̄dn |τii−1 > Φ̄dn |τi .
τi−1
τ

That is, decreasing Φ̄dn |τii+1 without increasing Φ̄dn |ττii−1 is impossible. This result also holds for
any interval [τi−1 , τi+1 ] (i = 1, · · · , Nτ − 1). Hence, we conclude that our proposed S-MGD
based method Pareto-optimally solves P2 . This completes the proof.
Remark 1. It is noteworthy that the proposed S-MGD method has polynomial computational
complexity, bounded by


3
3
O T ζτ,min
/µ ≤ ΘS ≤ O T 2 ζτ,max
/µ2 .

(48)

ζτ,min and ζτ,max are the minimal and maximal number of UAV-BSs in considered duration.
Proof. The number of time-slots during (0, T ) is T /µ. Accordingly, the number of greedy
searching shown in (46) is bounded by T /µ and T (T /µ + 1)/(2µ). Note that in each of these
greedy searchings, the optimal trajectory needs to be calculated by Hungarian method with a


τ
τi+1 3
complexity of O (ζτ |τi ) [25]. Since ζτ,min ≤ ζτ |τii+1 ≤ ζτ,max , (48) can be derived.
Remark 2. The sequential decision process shown in (46) is based on the global knowledge of
ground user density in (0, T ), which is inferred by the pattern formation systems.
Following Remark 2, we analyze the effects of density pattern accuracy on the performance
of proposed S-MGD method, as shown in following section.

20

V. L EARNING

THE

D ENSITY PATTERN

OF

G ROUND U SERS

This section shall focus on predicting the user density by machine learning techniques. In
particular, we first explicitly analyze the effect of pattern formation accuracy on the increment
of UAV-RF. Then, we characterize the feature of subregions that sufficient sampling is required
for reducing UAV-RF, based on the relation between sampling number and pattern accuracy.
A. Effects of Pattern Formation Accuracy
Recall that the minimum static UAV-RF is achieved at Rβ∗ , based on specific density of ground
users (conf. (29)). An inaccurate predicted density of users, denoted as λ̂β , would increase static
UAV-RF and reduce the life-time of batteries. Let the static UAV-RF corresponding to λ̂β be



Sβ Pcu
C/W
∗
∗2
+ λβ 2
− 1 Ptr,1 (hβ,1 )R̂β ,
Φ̂st,β (t) =
(49)
πEb R̂β∗2
where R̂β∗ is the generated coverage radius with λ̂β , i.e.
s
Pcu
R̂β∗ = 4
.
λ̂β (2C/W − 1) Ptr,1 (h∗β,1 )

(50)

The following theorem characterizes the average increased UAV-RF caused by inaccurate prediction on λβ .
Theorem 3. The expected value of increased static UAV-RF is
h
i
∆Φst,β (t) , E Φ̂st,β (t) − Φst,β (t) = Λβ (t)ξβ (t).

Λβ (t) is the eigenvalue of the β-th subregion at time t, given by
Sβ q
Λβ (t) =
Pcu (2C/W − 1) Ptr,1 (h∗β,1 )/λ3β .
4πEb

(51)

(52)

ξβ (t) , var(λ̂β (t)) + bias2 (λβ (t)) is the generalization error with quadratic loss function, which
is composed by the variance and the square of bias of predicted ground user density.
Proof. Since Rβ∗ minimizes Φst,β (t), ∂Φst,β (t)/∂Rβ |R∗β = 0 can be derived. Using Taylor series
expansion at Rβ∗ , Φ̂st,β (t) can be rewritten as
Φ̂st,β (t) =

Φ∗st,β (t)

1 ∂ 2 Φst,β (t)
+
2 ∂ 2 Rβ

Similarly, R̂β∗ (t) can be written as
R̂β∗ (t)

=

Rβ∗

∂Rβ
+
∂λβ

λβ



R∗β


h
i2
3 
∗
∗
∗
∗
R̂β − Rβ + O R̂β − Rβ
.

λ̂β − λβ




2 
+O
λ̂β − λβ
.

(53)

(54)

21

Substituting

∂ 2 Φst,β (t)
∂ 2 Rβ (t)

and

∂Rβ (t)
∂λβ (t)

into (53) and (54), one has

h
i
2 
E Φ̂st,β (t) − Φst,β (t) = Λβ (t)E λ̂β − λβ
,

where Λβ (t) is given by (52) and E
ξβ (t) , E



λ̂β − λβ



2 

λ̂β − λβ

=E

2 

h

(55)

is the generalization error defined as [26]


i2 
λ̂β − E(λ̂β ) + E(λ̂β ) − λβ


h
i2
= var(λ̂β ) + E(λ̂β ) − λβ

(56)

= var(λ̂β ) + bias2 (λβ ).

E(·) is the expected value of (·). Substituting (56) into (55), (51) can be proved.
Following Theorem 3, one can see that the generalization error of pattern formation systems
proportionally contribute to the expected value of increased static UAV-RF. Also, observing the
eigenvalue Λβ (t) shown in (52), we conclude that in subregions with small λβ , ∆Φst,β (t) is more
sensitive to generalization error. That is, when the time-varying density of ground users in one
specific subregion is low, the generalization error should be strictly restrained.
B. Bounds on Generalization Error and Sampling Number
Discussing the upper-bound of generalization error and the corresponding lower-bound of
sampling number under specific constraints is one of the core problems in the literature of
machine learning, since a large sample size leads to a lot of work [26], [27]. Therefore, we
consider minimizing the total sampling number of system under the constraint of maximal static
UAV-RF increment of considered area.
According to Vapnik-Chervonenkis (VC) theorem [20], the generalization error convergences
with respect to increasing sample size, i.e.
ξβ (t) ≤ ξˆβ (t) + ε(d, Ns,β (t), δ)
holds true with probability at least 1 − δ, where
s
ε(d, Ns,β (t), δ) =



1
log d + log
,
2Ns,β (t)
δ
1

(57)

(58)

ξˆβ (t) is the training error, d is the volume of hypothesis space and δ ∈ (0, 1). Ns,β (t) is the
sampling number in the β-th subregion at time t. Considering the fact that training error is

22

determined by the capacity of learning algorithm [26], we let ξmax be the maximal training error
and (57) can be rewritten as
ξβ (t) ≤ ξmax + ε (d, Nβ,s (t), δ) .

(59)

Denote the maximal tolerable static UAV-RF increment of considered area as ∆Φst,max (t).
Then, the optimization problem on minimizing the sampling number of considered area can be
expressed as follows.
P3 :

min
{ξβ (t)}

κ
X

Nβ,s (t)

(60)

β=1

s.t. ξβ (t) ≥ 0,
κ
X

(60a)

∆Φst,β (t) ≤ ∆Φst,max (t).

(60b)

β=1

(60a) denotes the non-negativity of generalization error, and (60b) shows the upper-bound of
the average increment of static UAV-RF in considered area at time t. By solving P3 , the
minimal sampling number and corresponding maximal generalization error of each subregion in
considered area are presented as Proposition 1.
Proposition 1. With maximal tolerable increment of static UAV-RF in considered area ∆Φst,max (t),
the minimal sampling number in the β-th subregion at t is
Nβ,s (t) = ω 2 Λ2β (t) (log d − log δ) /8,

(61)

where
ω=

2κ
P
.
∆Φst,max (t) − ξmax κβ=1 Λβ (t)

(62)

Correspondingly, the maximum generalization error of pattern formation system is
ξβ (t) ≤ ξmax +

2
.
ωΛβ (t)

Proof. The Lagrangian function of P3 is
!
!
κ
κ
κ
X
X
X
∆Φst,β (t) − ∆Φst,max (t) .
Nβ,s (t) − ω
Nβ,s (t), ω =
L
β=1

Let ∂L

P
κ

β=1

(63)

(64)

β=1


β=1 Nβ,s (t), ω /∂Nβ,s (t) = 0, we have
q
2 = ωΛβ (t) (log d − log δ) / (2Nβ,s (t)),

(65)

23

TABLE II
C OMMON PARAMETERS USED IN NUMERICAL ANALYSES

urban (a, b, η0 , η1 )

dense urban (a, b, η0 , η1 )

suburban (a, b, η0 , η1 )

bandwidth W

noise density N0

(9.61, 0.16, 1, 20)

(12.08, 0.11, 1.6, 23)

(4.88, 0.43, 0.1, 21)

10 KHz

5×10−15 W/Hz

carrier frequency f

BS coverage area Sbs

velocity vh , va , vd

data rate C

considered area size

2.4 GHz

104 m2

1 m/s

10 Kbps

1000 m × 1000 m

which is equivalent to (61). Since ω is the positive value that subjects to

Pκ

β=1 ∆Φst,β (t)

=

∆Φst,max (t), (62) can be derived. Also, one can obtain (63) by substituting (61) into (59).
With Proposition 1, we can conclude that in our considered area, subregions with high
eigenvalue need large sampling number. Further, let κ = 1, Proposition 1 implies that
Ns,1 (t) =

Λ2β (t) (log d − log δ)
2 (∆Φst,max (t) − ξmax Λβ (t))2

,

(66)

and
ξ1 (t) ≤ ∆Φst,max (t)/Λβ (t).

(67)

This agrees with our previous result that ∆Φst,β (t) in subregions with small λβ is more sensitive
to generalization error, as shown in Theorem 3.
VI. N UMERICAL R ESULTS
In this section, we present some numerical results to validate our theoretical founds. Besides,
more insights on the effectiveness of proposed optimal placement strategy and the corresponding
S-MGD based placement updating method are provided. The common parameters are listed in
Table II [6]. In our simulations, without specification, the simulated communication environment
is Urban. The considered area is a square, where two rectangle-shaped subregions with equal area
are investigated. A RSC is deployed at the center of considered area. Without loss of generality,
let Sβ /(πEb ) = 1 m2 /J, β = 1, 2 and Ph , Pa , Pd = Pm , where Pm is the mobility power.
A. Optimal Placement Strategy in Single Time-slot
Fig. 5 depicts the hovering altitude versus coverage radius in various environments when the
transmit power is fixed. In dense urban, urban and suburban environments, the red-solid lines
correspond to transmit power 0.05 W, 0.5 W and 1.5 W. Similarly, the blue-dash lines correspond

24

1000

R h * ,1

Hovering altitude [m]

800

LP
HP

Dense
Urban

600

h * ,1

h * ,1

=1.16

=1.86

400

h * ,1=0.43

200
Suburban

Urban

0
0

200

400

600

800

1000

Coverage radius [m]

Fig. 5. The hovering altitudes versus Rβ in various environments. The optimal hovering altitude is marked by stars.
10 6

3.5

800

P cu=0.5 W

3

P cu=5 W

700

P cu=50 W

600

(t)

2

1.5

=1 m -2,

R 2 =184.05 m

500

2

y [m]

Static UAV-RF

2.5

*
st,

400
3

=5 m -2,

R 3 =123.08 m

300
1
200
1

0.5

=0.1 m -2, R 1 =327.3 m

100

0

0
0

200

400

600

800

1000

1200

1400

Coverage radius [m]

Fig. 6. The static UAV-RF versus coverage radius with various
2

on-board circuit power. λβ =0.1 /m .

0

200

400

600

800

1000

x [m]

Fig. 7. An example of the coverage system under different
densities of ground users. Pcu =0.5 W.

to 0.15 W, 1.5 W and 5 W. The simulated density of ground users is 0.1 /m2 , and the optimal
hovering altitude is marked by stars. Observing the blue-dash line in suburban, one can find that
when the transmit power is fixed, the coverage radius achieves its maximum 810 m at hβ = 350
m. In other words, when the coverage radius is 810 m, hβ = 350 m is the optimal hovering
altitude that minimizes the transmit power. The solid black lines depict the optimal hovering
altitude with respect to Rβ . It can be seen that the optimal hovering altitudes is proportional
to the desired coverage radius, which verifies our theoretical results shown in Lemma 1. In
addition, it can be observed that in environment with high-rise buildings, the optimal hovering
altitude is also high. This is because high hovering altitude can reduce shade effects from high-

25

rise buildings, and hence the LOS probability is increased. In this way, the average path loss is
decreased and transmit power is saved.
The static UAV-RF versus coverage radius with various on-board circuit power is depicted in
Fig. 6. The density of ground users is 0.1 /m2 , and the optimal coverage radii that minimize
Φst,β (t) are marked by stars. When Pcu =0.5 W, 5 W and 50 W, the simulated Rβ∗ =327.3 m, 582
m and 1035 m. As expected, high on-board circuit power leads to high static UAV-RF, which
has been shown in (31). This indicates that restricting the on-board circuit power of UAV-BSs
can effectively decrease the static UAV-RF. Also, it can be observed that the optimal coverage
radius increases with respect to on-board circuit power. This is because when Pcu is high, large
Rβ (t) can decrease the number of UAV-BSs, resulting in the reduction of the total consumed
on-board circuit power of network.
Fig. 7 shows an example of coverage system under different densities of ground users, when
Pcu =0.5 W. The related user densities in three subregions are 0.1 /m2 , 1 /m2 and 5 /m2 ,
respectively. The corresponding theoretical coverage radii given by (28) are 327.3 m, 184.05 m
and 123.08 m, which agree with the simulated results. To the hovering altitudes, they can be
obtained with Lemma 1. It can be observed that in cases with dense users, the optimal coverage
is small. This is because small coverage radius can reduce the transmit power increased by high
user density. Comparing the simulation results in considered subregions, we can conclude that
our proposed optimal placement strategy can efficiently adjust to varying user densities, while
minimizing the static UAV-RF.
B. Optimal Placement Updating Method in Considered Duration
To demonstrate the effectiveness of our proposed S-MGD method, we illustrate two robust
placement updating methods for comparison. The first one is the ’Lazy’ method, where UAV-BSs
hold their placement in our considered duration. In contrast, The second one is the ’Diligent’
method, where UAV-BSs update their positions at the beginning of each time-slot. In this way,
the static UAV-RF of ’Diligent’ method, denoted as ΦD
st,β (t), is always minimized. When take
the energy cost on mobility into consideration, ΦD
st,β (t) can be denoted as the lower-bound of
Φdn,β (t). Here, we take ΦD
st,β (t) as the base line to evaluate the effectiveness of our proposed
method. In this subsection, Pcu = 0.5 W and T =24 h.
Fig. 8 depicts the updated coverage radii with our proposed S-MGD scheme, which reflects
the placement of UAV-BSs. Clearly, the number of update epochs depends on the mobility

26

power. Compare the cases with different mobility powers, one can observe that in cases with
low mobility power, the coverage radii of UAV-BSs updates frequently. By contrast, in cases
with high mobility power, UAV-BSs tend to keep their coverage radii for several time-slots.
Specifically, when Pcu =50 W, the coverage radii of UAV-BSs remain unchanged during [0, T ].
This is because updating coverage radii may cost more energy than it saves.

Coverage radius [m]

200
150

L
100
250

Base line
P m=0.05 W

D

P m=1.5 W

200

P m=50 W

150
0

8

16

24

32

40

48

56

64

72

Time [h]

Fig. 8. The coverage radii of UAV-BSs in cases with different mobility powers. Pcu = 0.5 W and T =24 h.

In addition, Fig. 8 shows that the placement updating frequency relies on the varying density
of ground users. For example, in the case with Pm =1.5 W, ’L’ and ’D’ labels the lazy and diligent
part of our proposed S-MGD scheme. It can be seen that in lazy part, the placement of UAV-BSs
remain unchanged. By contrast, in diligent part, the placement is updated in each time-slot. This
is intuitive, because in lazy part, the optimal coverage radii in both considered subregions return
to the initial values at the end of this part. Therefore, there is no need to update the placement.
Also, in diligent part, the optimal coverage radii vary sharply. Thus, spending energy on updating
the placement can save more energy cost on communication and on-board circuit. Simulation
results in Fig. 8 demonstrate that our proposed S-MGD based placement updating method can
efficiently adjust to the varying user density and mobility power.
Fig. 9 compares the dynamic UAV-RF of Diligent, Lazy and our proposed S-MGD based
placement updating method. The simulated mobility power is 0.05 W, 1.5 W and 50 W. Fig. 9
shows that S-MGD always provides the most energy efficiency updating strategy for arbitrary
mobility powers. Also, as expected, in cases with high and low mobility power, the proposed
method degrades to ’Lazy’ and ’Diligent’ method, respectively. Observe the labeled average
dynamic UAV-RF, we can see that compared with the worst updating method, our method can
reduce the average dynamic UAV-RF by 7%∼96%.

27

Dynamic UAV-RF

2

1.5
40
1
0.6

Base line
S-MGD
Diligent
Lazy

0.5

0.5
0.4
24

0
0

8

16

24

32

32

40

48

56

64

72

t/h

Fig. 9. The UAV-RF in cases with different mobility powers. Pcu = 0.5 W and T =24 h.

In addition, since the ’Lazy’ method holds the placement of UAV-BSs during [0, T ], the
average dynamic UAV-RF remains unchanged for arbitrary mobility powers. Compared with the
base line, the ’Lazy’ updating method can provide satisfactory energy efficiency. However, Fig.
10 shows that the performance of ’Lazy’ method is unstable and may vary with the initial time of
our considered duration. By contrast, our method can stably provide near-optimal performance
in arbitrary cases, regardless of the initial time of considered duration and the mobility powers.
C. Insights on Learning Density Patterns
Fig. 11 compares the simulated and theoretical ∆Φst,β (t) with Monte-Carlo method, where the
repetition number is 106 . For the simplicity of comparison, the depicted ∆Φst,β (t) is normalized
by Φst,β (t). Obviously, our theoretical results shown in Theorem 3 are sufficiently explicit in
most of the simulated cases. However, the difference between our derived theoretical results and
the simulated ∆Φst,β (t) grows with respect to generalization error. This is because Proposition
1 is derived with Taylor series expansion, which lost accuracy for large λ̂β − λβ . In addition,
compare the results shown in cases with λβ =3 m−2 and 10 m−2 , we can see that our theoretical
results in Theorem 3 are more accurate when user density is high. This is because as previous
illustrated, ∆Φst,β (t) in subregions with low user density are more sensitive to generalization
error.
The lower-bound of sampling numbers in considered area and subregion versus user density
λβ are depicted in Fig. 12(a) and Fig. 12(b), respectively. In simulations, ∆Φst,max (t)=10, and
without loss of generality, let ξmax =0. For our considered subregions, the simulated density of
users is λβ and 1 − λβ , where λβ ∈ [0, 1/2]. Also, some sampling numbers satisfying (60b)
are labeled with Monte-Carlo method. It can be observed that the sampling numbers are lower-

28

2.2
1.5
2

S-MGD
Diligent
Lazy

1.4

S-MGD
Diligent
Lazy

1.8

1.3
1.6
1.2

1.4

1.1

1.2

1

1
0

6

12

18

24

Time [h]

0

6

12

18

24

Time [h]

Fig. 10. The average dynamic UAV-RF of UAV-BSs versus different initial time.
10 -3
Monte Carlo
Theoretical

Average increment of static UAV-RF

Average increment of static UAV-RF

0.02

0.015

0.01

0.005

0

Monte Carlo
Theoretical

7
6
5
4
3
2
1
0

0

1

2

0
-2

Generalization error [m ]

2

4

Generalization error [m -2 ]

Fig. 11. The simulated and theoretical expected increase of
6

static UAV-RF versus user density. Repetition number 10 .

Fig. 12. The lower-bound of sampling numbers in considered
area (a) and the first subregion (b) versus user density λβ .

bounded by (61), which verifies Proposition 1. Also, Fig. 12(b) shows that cases with low user
density require large sample sets. This is consistent with the intuition that in cases with lower
user densities, ∆Φst,β (t) is more sensitive to generalization error. This manifests that cases with
low user densities may effect the energy performance of network greatly and need to be taken
into consideration carefully.
VII. C ONCLUSION
This paper focused on the placement of UAV-BSs for the downlink system, where the powers
relevant to signal transmitting, on-board circuit and the potential mobility of UAVs as well as the
density of ground users are taken into account. To track the instable and non-ergodic time-varying

29

nature of user density, we first proposed a pattern formation based framework in a machine
learning manner. Then, to minimize the UAV-RF of network, the optimal placement of UAVBSs was analyzed. In the static case with one time-slot, we proved that the optimal placement is
achieved when the transmit power of UAV-BSs equals their on-board circuit power. In addition,
in the dynamic case with a multiple time-slot duration, we showed that the optimal placement updating problem is an nonlinear dynamic programming coupled with an inherent integer
linear programming. Considering the NP-hardness of original problem, we proposed a S-MGD
based Pareto-optimal placement updating method with a polynomial time complexity. Simulation
results showed that it can stably provide near-optimal performance in terms of minimum UAVRF. Finally, we proved that the generalization errors of learned density patterns proportionally
contribute to the increment of UAV-RF in each subregion. Simulation and theoretical results
demonstrated that large sample sets are needed in regions with large area, high-rise buildings or
low user density.
Further, it is seen from Remark 1 that the computational complexity of proposed S-MGD
method is greatly affected by the number of UAV-BSs. In previous mentioned regions, the computational complexity can be extremely high. Therefore, to reduce the complexity, multiple RSCs
and proper partition to considered region are needed. In addition, the effect of pattern formation
accuracy on the dynamic UAV-RF should be discussed. These issues shall be considered in future
works.
R EFERENCES
[1] L. Gupta, R. Jain, and G. Vaszkun, “Survey of important issues in UAV communication networks,” IEEE Commun. Surveys
Tuts., vol. 18, no. 2, pp. 1123–1152, 2016.
[2] M. Mozaffari, W. Saad, M. Bennis, and M. Debbah, “Efficient deployment of multiple unmanned aerial vehicles for optimal
wireless coverage,” IEEE Commun. Lett., vol. 20, no. 8, pp. 1647–1650, 2016.
[3] E. Kalantari, H. Yanikomeroglu, and A. Yongacoglu, “On the number and 3D placement of drone base stations in wireless
cellular networks,” in Proc. IEEE Vehicular Technology Conf., 2016.
[4] M. Mozaffari, W. Saad, M. Bennis, and M. Debbah, “Optimal transport theory for power-efficient deployment of unmanned
aerial vehicles,” in Proc. IEEE ICC, 2016, pp. 1290–1295.
[5] A. Al-Hourani, S. Kandeepan, and S. Lardner, “Optimal lap altitude for maximum coverage,” IEEE Wireless Commun.
Lett., vol. 3, no. 6, pp. 569–572, 2014.
[6] R. I. Bor-Yaliniz, A. El-Keyi, and H. Yanikomeroglu, “Efficient 3-d placement of an aerial base station in next generation
cellular networks,” in Proc. IEEE ICC, 2016, pp. 985–989.
[7] J. Lyu, Y. Zeng, R. Zhang, and T. J. Lim, “Placement optimization of UAV-mounted mobile base stations,” IEEE Commun.
Lett., vol. 21, no. 3, pp. 604–607, 2017.

30

[8] M. Alzenad, A. El-Keyi, F. Lagum, and H. Yanikomeroglu, “3D placement of an unmanned aerial vehicle base station
(UAV-BS) for energy-efficient maximal coverage,” IEEE Wireless Commun. Lett., vol. 6, no. 4, pp. 434–437, 2017.
[9] C. Di Franco and G. Buttazzo, “Energy-aware coverage path planning of UAVs,” in Proc. IEEE ICARSC, 2015, pp.
111–117.
[10] J. Lu, S. Wan, X. Chen, and P. Fan, “Energy-efficient 3D UAV-BS placement versus mobile users’ density and circuit
power,” in Proc. IEEE GlobeCom, 2017.
[11] M. Mozaffari, W. Saad, M. Bennis, and M. Debbah, “Mobile unmanned aerial vehicles (UAVs) for energy-efficient Internet
of things communications,” arXiv preprint arXiv:1703.05401, 2017.
[12] M. Asllani, J. D. Challenger, F. S. Pavone, L. Sacconi, and D. Fanelli, “The theory of pattern formation on directed
networks,” Nature Commun., vol. 5, p. 4517, 2014.
[13] V. Sharma, U. Mukherji, V. Joseph, and S. Gupta, “Optimal energy management policies for energy harvesting sensor
nodes,” IEEE Trans. Wireless Commun., vol. 9, no. 4, 2010.
[14] E. Baştuğ, M. Bennis, E. Zeydan, M. A. Kader, I. A. Karatepe, A. S. Er, and M. Debbah, “Big data meets telcos: A
proactive caching perspective,” Journal of Commun. Netw., vol. 17, no. 6, pp. 549–557, 2015.
[15] M. Bkassiny, Y. Li, and S. K. Jayaweera, “A survey on machine-learning techniques in cognitive radios,” IEEE Commun.
Surveys Tuts., vol. 15, no. 3, pp. 1136–1159, 2013.
[16] M. A. Alsheikh, S. Lin, D. Niyato, and H. P. Tan, “Machine learning in wireless sensor networks: Algorithms, strategies,
and applications,” IEEE Commun. Surveys Tuts., vol. 16, no. 4, pp. 1996–2018, 2014.
[17] M. Kulin, C. Fortuna, E. De Poorter, D. Deschrijver, and I. Moerman, “Data-driven design of intelligent wireless networks:
An overview and tutorial,” Sensors, vol. 16, no. 6, 2016.
[18] G. M. Dias, B. Bellalta, and S. Oechsner, “A survey about prediction-based data reduction in wireless sensor networks,”
ACM Comput. Surveys, vol. 49, no. 58, 2016.
[19] K.-S. Lee, S.-R. Lee, Y. Kim, and C.-G. Lee, “Deep learning–based real-time query processing for wireless sensor network,”
International J. Distrib. Sensor Netw., vol. 13, no. 5, 2017.
[20] V. Vapnik, The nature of statistical learning theory.

Springer science & business media, 2013.

[21] H. Wang, F. Xu, Y. Li, P. Zhang, and D. Jin, “Understanding mobile traffic patterns of large scale cellular towers in urban
environment,” IEEE/ACM Trans. Netw., 2016.
[22] M. Mozaffari, W. Saad, M. Bennis, and M. Debbah, “Unmanned aerial vehicle with underlaid device-to-device communications: Performance and tradeoffs,” IEEE Trans. Wireless Commun., vol. 15, no. 6, pp. 3949–3963, 2016.
[23] A. Srinivas, G. Zussman, and E. Modiano, “Construction and maintenance of wireless mobile backbone networks,”
IEEE/ACM Trans. Netw., vol. 17, no. 1, pp. 239–252, 2009.
[24] H. A. Fayed and A. F. Atiya, “A mixed breadth-depth first strategy for the branch and bound tree of euclidean k-center
problems,” Comput. Opti. Appl., vol. 54, no. 3, pp. 675–703, 2013.
[25] H. W. Kuhn, “The hungarian method for the assignment problem,” Naval research logistics quarterly, vol. 2, no. 1-2, pp.
83–97, 1955.
[26] I. Goodfellow, Y. Bengio, and A. Courville, Deep learning.

MIT press, 2016.

[27] R. Dobrusin, “A statistical problem arising in the theory of detection of signals in the presence of noise in a multi-channel
system and leading to stable distribution laws,” Theory Prob. Its Appl., vol. 3, no. 2, pp. 161–173, 1958.

