A QoS-Oriented Trajectory Optimization in Swarming Unmanned-Aerial-Vehicles Communications

TL;DR

This paper introduces a new trajectory and resource allocation optimization method for UAV swarms to enhance communication quality, balancing fairness among ground stations through an iterative convex approximation approach.

Contribution

It presents a novel joint optimization framework for UAV trajectories, power, and scheduling to meet QoS requirements, using successive convex approximation and block coordinate techniques.

Findings

Improved fairness in UAV communication networks.

Effective iterative algorithm for joint trajectory and resource optimization.

Enhanced QoS adherence in UAV-based systems.

Abstract

This letter aims to present a novel approach for unmanned aerial vehicles (UAV)' path planning with respect to certain quality of service requirements. More specifically, we study the max-min fairness problem in an air-to-ground communication system where multiple UAVs and multiple ground stations exist. We jointly optimize the UAVs trajectories and power allocation as well as the user scheduling. To this end, we propose an effective iterative algorithm that relies on the successive convex approximation and the block coordinate decent techniques.

Tables

Table 1. Table I : Convergence Behavior

$r$	1	2	3	4	-	8	9	10	11
$\Delta {\mu }^{lb}$	125	16.5	8	2.3	-	0.3	0.2	0.2	0.1

Table 2. Table II : Performance achieved by the different ground stations

GT No.	1	2	3	4	5	6
Data rate	209.1	208.6	209.4	210.1	208.8	208.4
(Bps/Hz)
Connection time	34	31	31	32	36	31
(sec)

Eq.	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
$b$	$M\left(N-1\right)$	$M\left(N-1\right)$	$MN$	$MN$	$K$	$KMN$	$M$	$MN$	$MN$	$\frac{M\left(M-1\right)N}{2}$

Equations

$${\gamma }_{k,m}(n)=\frac{p_m(n)h_{k,m}(n)}{{\sum }_{\begin{array}{c}j=1\\ j\ne m\end{array}}^M{p}_j(n)h_{k,j}(n)+{\sigma }_0^2},$$

$$R_k=\frac{1}{N}\sum _{n=0}^N\sum _{m=1}^M{\alpha }_{k,m}(n){\log }_2(1+{\gamma }_{k,m}(n)).$$

$$E_m$$

$${\mathbf{v}}_m(n+1)$$

$${\mathbf{q}}_m(n+1)$$

$$\Vert {\mathbf{v}}_m(n)\Vert$$

$$\Vert {\mathbf{a}}_m(n)\Vert$$

$$\Vert {\mathbf{v}}_m(n)\Vert$$

$$\Vert {\mathbf{q}}_m(n)-{\mathbf{q}}_j(n)\Vert$$

$$\sum _{k=1}^K{\alpha }_{k,m}(n)\le 1,\text{indent}\sum _{m=1}^M{\alpha }_{k,m}(n)\le 1,\text{indent}\forall k,m,n.$$

$$\underset{\mu ,\mathcal{Q},\mathcal{V},\mathcal{A},\mathcal{S},\mathcal{P}}{\max }$$

$$(\text{P1})$$

$$\underset{\mu ,\mathcal{Q},\mathcal{V},\mathcal{A},\mathcal{S},\mathcal{P}}{\max }$$

$$(\text{P2})$$

$$(\text{P2.1})$$

$$(\text{P2.2})$$

$$R_k=\frac{1}{N}\sum _{n=0}^N\sum _{m=1}^M{\alpha }_{k,m}(n){\log }_2\left(1+\frac{B_{k,m}(n)}{{\sum }_{\begin{array}{c}j=1\\ j\ne m\end{array}}^M{B}_{k,j}(n)+{\sigma }_0^2}\right).$$

$$0\le B_{k,m}(n)\le \frac{{\beta }_0{P}_{max}}{H^2+\Vert {\mathbf{q}}_m(n)-{\mathbf{w}}_k{\Vert }^2}\text{indent}\forall k,m,n,$$

$$(\text{P3})$$

$$\underset{{\mu }^{lb},\mathcal{B},\mathcal{Q},\mathcal{V},\mathcal{A},\Lambda }{\max }$$

$$D_{k,m}(n)$$

$$F_{k,m}(n)$$

$${\overline{R}}_{k,m}^{upper}(n)$$

$$\underset{{\mu }^{lb},\mathcal{B},\mathcal{Q},\mathcal{V},\mathcal{A},\Lambda }{\max }$$

$$\Vert {\mathbf{v}}_m(n){\Vert }^2\ge \Vert {\mathbf{v}}_m^r(n){\Vert }^2+2{\left({\mathbf{v}}_m^r(n)\right)}^T({\mathbf{v}}_m(n)-{\mathbf{v}}_m^r(n)).$$

$${\lambda }_m(n{)}^2\le \Vert {\mathbf{v}}_m^r(n){\Vert }^2+2{\left({\mathbf{v}}_m^r(n)\right)}^T({\mathbf{v}}_m(n)-{\mathbf{v}}_m^r(n)),\forall m,n$$

$$R_{k,m}(n)$$

$$A_{k,j}(n)=$$

$$C_{k,j}(n)=$$

$$\frac{1}{N}\sum _{n=0}^N\sum _{m=1}^M\left({\alpha }_{k,m}(n){\log }_2(\sum _{\begin{array}{c}j=1\end{array}}^M{B}_{k,j}(n)+{\sigma }_0^2)-{\overline{R}}_{k,m}^{upper}(n)\right)\ge \mu$$

$$B_{k,m,max}(n)\ge P_{max}{\beta }_0\left(-\frac{\Vert {\mathbf{q}}_m(n)-{\mathbf{w}}_k{\Vert }^2}{H^4}+D_{k,m}(n)({\mathbf{q}}_m(n)-{\mathbf{w}}_k{)}^T({\mathbf{q}}_m^r(n)-{\mathbf{w}}_k)+F_{k,m}(n)\right),$$

$$D_{k,m}(n)$$

$$F_{k,m}(n)$$

$$0\le B_{k,m}(n)$$

$$\mu ({\mathcal{S}}^r,{\mathcal{Q}}^r,{\mathcal{V}}^r,{\mathcal{B}}^r)\le \mu ({\mathcal{S}}^{r+1};{\mathcal{Q}}^r,{\mathcal{V}}^r,{\mathcal{B}}^r)$$

Full text

I Introduction

Recently, unmanned aerial vehicles (UAVs) have drawn significant research interests [1, 2]. This is mainly due to the significant decrease in the production costs of these systems [3], added to their flexible deployment and high maneuverability. Traditionally, the use of UAVs was only dedicated to military applications in the late 90’s [4]. However, research and development of UAVs have enabled the emergence of several new applications. As a matter of fact, the drones have started to be used in cargo delivery, traffic control, weather forecasting, farming [5], etc.

Motivated by this context, earlier research started by studying single UAV systems. To begin with, they studied the problem of path planning that consists into scheduling the areas to visit in a given time period. For instance, the authors in [6] presented a multi-colony algorithm that optimizes the UAV path planning with obstacle avoidance.

After that, recent research has investigated the use of multi-UAV systems. In fact, it has been proved that, compared to a single UAV system, the use of a fleet of UAVs allows to perform more complex and challenging missions with higher speed and especially with increased performance and efficiency. However, with the use of a multi-UAV system, we are facing a new kind of problems such as interference, energy efficiency, path planning and user scheduling. For instance, the work in [7] has treated the path planning problem for a multi-UAV system using genetic algorithms and Bezier curves. Furthermore, a game theoretic approach has been presented in [8] to reduce the interference between the UAVs. However, there is only the work in [9] that considered working on all these problems combined. However, in that work, every single one of these problem was optimized separately. In the same context, we will conduct different optimization approach to get better fairness among served users.

In this letter, we study a multi-UAV aided wireless system where a group of UAVs are used to provide down-link communication to a group of ground stations (GT) in a two dimensional area given certain quality of service requirements in communication. The objective here is to maximize the minimum average rate among ground stations in order to achieve fair performance among them during a given finite period. To this extent, we consider a joint optimization problem in which we have to provide an efficient flight trajectory solution for every UAV with specifying the optimal UAV-user association and power allocation at each time. To that end, we propose an iterative algorithm based on the successive convex approximation (SCA) and the block coordinate descent (BCD) techniques. The proposed algorithm and the resulting fairness will be verified by simulation results.

II System model and problem formulation

We consider a wireless communication system where $M$ UAVs are employed to communicate with $K$ ground stations. Furthermore, we assume that every ground station $k\;(k=1,..,K)$ has a fixed position $\mathbf{w}_{k}=[x_{k},\;y_{k}]$ on the ground, while all the UAVs $(m=1,..,M)$ fly on the two dimensional coordinate $\mathbf{q}_{m}(n)=[x_{m}(n),\;y_{m}(n)]$ at a fixed altitude $H$111We suppose here that the trajectory optimization problem is intractable over continuous time. For ease of exposition, we use the discrete linear state-space approximation. This method consists of discretizing the time horizon $T$ into $N$ time slots with a step size $\delta t$. So, we have $t=n\delta t$ with $n=0,1,\ldots,N$. Thus all the formulas in this letter will be defined as a function of the discrete time instant $n$.. Thus, the distance between an UAV $m$ and a ground station $k$ is expressed by $d_{k,m}(n)=\sqrt{H^{2}+\|\mathbf{q}_{m}(n)-\mathbf{w}_{k}\|^{2}}$. For ease of exposition, we assume that the high altitude of UAVs enables them to effectively establish line-of-sight (LoS) link and the Doppler effect caused by their mobility is well compensated at the receivers. Therefore, the time-varying channel between an UAV $m$ and a ground station $k$ follows the free-space path loss model, which can be expressed as $h_{k,m}(n)=\frac{\beta_{0}}{H^{2}+\|\mathbf{q}_{m}(n)-\mathbf{w}_{k}\|^{2}}$, where $\beta_{0}$ represents the channel gain at reference distance ($d=1$ m).

We also assume that all the UAVs operate in the same frequency band. Thus, the corresponding received signal to interference plus noise ratio (SINR) at user $k$ from UAV $m$ is computed as

[TABLE]

where $p_{m}(n)$ denotes the down-link transmit power of UAV $m$, $\sigma_{0}^{2}$ represents the variance of additive white Gaussian noise (AWGN) at the receiver and the term $\sum_{\begin{subarray}{c}j=1\\ j\neq m\end{subarray}}^{M}p_{j}(n)h_{k,j}(n)$ stands for the co-channel interference (CCI) received by ground station $k$ during its communication with the UAV $m$ at time $n$. We consider the power control of the signals transmitted by UAVs to ground stations at each time $n$ in order to manage co-channel interference to the other users. The down-link transmit power of an UAV $m$ at time instant $n$ is constrained on $0\leq p_{m}(n)\leq P_{max}$, where $P_{max}$ denotes the maximum allowable transmit power at the UAV.

In order to define the UAV-user association and scheduling, we use a binary variable $\alpha_{k,m}(n)$ that takes 1 during communication between an UAV $m$ and a ground station $k$ in time $n$. Otherwise, it is equal to 0. Then, we can define the achievable rate of user $k$ over the period $N$ is given by:

[TABLE]

Each UAV operates with a finite amount of energy due to the compact size requirement which is limiting the serving time and data rate on flight. As explained in [123], we assume that the communication energy consumption is much smaller than the propulsion energy consumption and we only consider the propulsion power as the source of energy consumption, which can be expressed by

[TABLE]

where $c_{1}$ and $c_{2}$ are constants, $g$ is the gravitational acceleration, $w$ is the mass of an UAV and $E_{max}$ is the maximum allowable propulsion energy during the period $N$. In (3) $\mathbf{v}_{m}(n)$ and $\mathbf{a}_{m}(n)$ are the velocity and acceleration of UAV $m$ at time instant $n$ which are approximated by the first order Taylor expansion as

[TABLE]

The velocity and acceleration speed are limited to

[TABLE]

where $V_{max}$ and $V_{min}$ are the maximum and minimum velocity during flight and $a_{max}$ is the maximum acceleration.

In order to avoid collusion, we add a constraint that sets a minimum distance between the UAVs which is defined as

[TABLE]

We suppose that, at each time slot, each UAV serves at most one ground station and that each ground station is served at most by one UAV. Thus, the constraints on association binary variables are given by

[TABLE]

III Max-Min Rate Problem Formulation

In order to better understand the problem, we adopt the following notations :

• •

$\mathcal{S}$= $\{\alpha_{k,m}(n)$   $\forall k,m,n\}$: the user-UAV association.

• •

$\mathcal{P}$= $\{p_{m}(n)$,   $\forall m,n\}$: the transmit power of the UAVs.

• •

$\mathcal{Q}$= $\{\mathbf{q}_{k,m}(n)$,   $\forall k,m,n\}$: the UAVs’ trajectories.

• •

$\mathcal{V}$= $\{\mathbf{v}_{k,m}(n)$,   $\forall k,m,n\}$: the UAVs’ velocities.

• •

$\mathcal{A}$= $\{\mathbf{a}_{k,m}(n)$,   $\forall k,m,n\}$: the UAVs’ accelerations.

As explained in the first section, we want to jointly optimize the down-link communication and the UAV-user scheduling as well as the trajectory, velocity and acceleration of each UAV in order to maximize the minimum average rate among all users. To this end, if we set $\mu=\min\limits_{k\in\mathcal{K}}R_{k}$, we end up with the following problem (P1) :

[TABLE]

Problem (P1) is challenging to solve due to the following reasons. To begin with, our problem includes integer constraint that is expressed in (13) due to the presence of the binary optimization variables of UAV-user scheduling and association. Furthermore, the constraints (3), (8), (9), and (B-B) are not convex with respect to either $\mathcal{Q,V,A,P}$. Finally, all the optimization variables are closely correlated. For instance, even a slight change of a given UAV trajectory can have an enormous impact on the trajectories of the other UAVs. So, we need to jointly optimize these variables in order to reach the better solution. All these constraints yield a mixed-integer non-convex problem and, in general, there is no standard method for solving this kind of problems efficiently.

IV Proposed Algorithm

Our problem is a mixed-integer non-convex problem that involves integer constraints. To tackle our problem of interest, we relax our problem by changing the binary variables of UAV-user scheduling into continuous variables. So, the resulting problem (P2) is:

[TABLE]

This relaxation implies that the objective value of the new problem (P2) will be considered as an upper bound of the problem (P1). Thus, from now on, we will optimize this new problem (P2) which is still non-convex due to the constraints (3), (8) and (B-B).

In order to do so, we use an iterative algorithm by leveraging the BCD method. Specifically, we divide our optimization variables into two blocks; the first block is for the user scheduling and association ($\mathcal{S}$) and the second one is represented by the UAVs trajectory and power control ($\mathcal{Q,V,A,P}$). So, these two blocks are alternately optimized in each iteration. As a result, we obtain two sub-problems of problem (P2).

In the first sub-problem, we only retain the constraints that are related to the user scheduling, while other variables ($\mathcal{Q,V,A,P}$) are fixed. Thus the first sub-problem (P2.1) can be expressed as:

[TABLE]

(P2.1) is a standard linear problem. So, it can be solved directly by using some known optimization tools.

In the second sub-problem (P2.2), we only retain the constraints that depend on the UAVs’ trajectory, velocity, acceleration and power control. Thus, (P2.2) is defined by

[TABLE]

The sub-problem (P2.2) is still non-convex due to the constraints (3), (8), (9) and (B-B). In order to solve it, we first couple the trajectory and the transmit power variables by introducing an auxiliary variable $B_{k,m}(n)$ which is defined by $B_{k,m}(n)=p_{m}(n)h_{k,m}(n)$. As a result, $R_{k}$ is rewritten as

[TABLE]

Additionally, the constraint (13) is transformed into

[TABLE]

which is still non-convex. Thus, if we set $\mathcal{B}$= [$B_{k,m}(n),$ $\forall k,m,n$], we have the new problem (P3) which can be formulated as

[TABLE]

In the following, we will try to solve the sub-problem (P3) using the successive convex optimization techniques. The main idea consists in finding another convex problem which represents a lower bound of (P3). Then, we successively maximize the lower bound of (P3) via optimizing the trajectory and the auxiliary variable.

Using the approximations in (Appendix A), we can find a convex problem (P4) for a given trajectory, velocity and the auxiliary variable ($\mathcal{Q}^{r},\mathcal{V}^{r},\mathcal{B}^{r}$) at iteration $r$ defined by:

[TABLE]

where $\Lambda=\{\lambda_{m}(n)$   ,$\forall m,n$} is a set of slack variables that we introduced to the problem, and the constants $D_{k,m}(n)$, $F_{k,m}(n)$ and $\overline{R}_{k,m}^{upper}(n)$ are given by:

[TABLE]

It is worthy noting that the feasible set of problem (P4) can be considered as a subset of the sub-problem (P2.2). Thus, by solving iteratively the problem (P4), we will be successively optimizing the lower bound of sub-problem (P2.2) in order to better approximate it. Thus, we will end up with an efficient approximate solution of the latter. Finally, we present an overall algorithm that synthesizes all the work done in this perspective:

  **Algorithm : (BCD) Algorithm for Problem (P2). **

Input : Locations of GTs, $V_{max},V_{min},A_{max},P_{max},$ $E_{max},\sigma_{0}^{2},\beta_{0},c_{1},c_{2}$ and $\epsilon$

Initialize $\mathcal{Q}^{0},\mathcal{V}^{0},\mathcal{B}^{0}$ and $r=0$ 2. 2.

Repeat

1. (a)

Solve problem (P2.1) for given $\mathcal{Q}^{r}$,$\mathcal{V}^{r}$,$\mathcal{B}^{r}$, and denote the optimal solution as $\mathcal{S}^{r+1}$. 2. (b)

Solve problem (P4) for given $\mathcal{S}^{r+1}$, and denote the optimal solutions as $\mathcal{Q}^{r+1}$,$\mathcal{V}^{r+1}$,$\mathcal{B}^{r+1}$. 3. (c)

Update $r=r+1$. 3. 3.

Until   The fractional increase of the objective value of problem (P2) is below a threshold $\epsilon$ . 4. 4.

Define   $p_{m}(n)=\displaystyle\sum_{k=1}^{K}\alpha_{k,m}(n)\frac{B_{k,m}(n)}{h_{k,m}(n)},\;\;\;\forall m,n$.

As explained in (Appendix B), our algorithm performs better at each iteration and is guaranteed to converge. Also, the complexity of our algorithm is of order $\Theta\left((KM)^{3}N+(KMN)^{2}(KMN+\frac{M^{2}N}{2})\right)$. But due to space limitation, the proof of convergence and complexity are shown in the companion report [10]. Also, it has been proved in [11] that the convergence complexity of our algorithm in terms of number of iterations is $\Theta(1/r)$ (i.e.$f\left(x^{r}\right)-f^{*}=\Theta\left(\frac{1}{r}\right)$, where $f^{*}$ is the optimal objective value and $f\left(x^{r}\right)$ is the objective value at iteration r). Thus, this indicates that our algorithm converges in a few iterations. We will prove that our algorithm performs very well by testing its effectiveness on some examples.

V Numerical results

We consider a system model composed of $M=2$ UAVs at a fixed altitude $H=100$ m and $K=6$ GTs on the ground that are arbitrarily distributed in an area of $500\times 500$ m2 for simulation purpose. For channel gain, we set $P_{max}=0.1$ Watt, $\beta_{0}=-60$ dB and $\sigma_{0}^{2}=-110$ dBm, respectively. For propulsion energy and flight dynamics, we set $V_{min}=1.5$ m/s, $V_{max}=50$ m/s, $a_{max}=5$ m/s2, $E_{max}=2\times 10^{5}J$, $c_{1}=9.26\times 10^{−4}$ and $c_{2}=2250$, respectively. Simulation period is $N=100$ with step size $\delta t=1$ sec.

Next, in order to use our algorithm, we need to provide an initial scheme for the different variables of our model. To that end, we start by providing a simple systematic circular trajectory design for all the UAVs using the $K$-mean clustering algorithm to find the radius and center (also called centroid) of each trajectory as explained in (Appendix C). We also assume that the initial speed of every UAV $m$ is set to be constant during all the flight which are respectively equal to $3$ and $4m/s$.

Regarding the user scheduling, at each time slot, every UAV $m$ will communicate with the closest ground station that belongs to its assigned cluster $m$ that results from the $K$-mean clustering algorithm. Also, the transmission power of all the UAVs is set to be at its maximum level $P_{max}$.

Thus, we end up with an efficient initialization scheme for our problem due to the following reasons. First, this method offers a feasible solution that is compliant with all the constraints of our problem. Second, by using this method, we ensure that all the UAVs will provide full coverage for the areas that contains the ground stations. As a result, all the ground stations will be served equally during the flight. Thus, this will help achieve fair performance between all the ground stations which represents our main goal from the beginning.

Finally, since the $K$-mean clustering algorithm divides the ground stations into clusters that are generally far away from each other, we can guarantee that the UAVs are sufficiently separated during their communications with the ground stations. As a result, this tends to minimize the co-channel interference on the ground stations and thus, the achievable data rate at the ground stations will be at its maximum level.

We note that, in the initial solution that we provided, the UAV 2 communicates only with two ground stations which are represented by the red dots in Figure 1(a). After applying our proposed algorithm, the UAV 2 extended its trajectory in order to serve a third ground station that was only served by the first UAV in the initial solution as shown in Figure 1(a).

We also note that every UAV tries to position itself near the ground station when there is a communication between them. This is explained by the fact that this helps them achieve better communication channel. For instance, all the UAVs tend to have a high velocity when they are moving from a ground station to another and tend to have a minimum velocity when they are hovering next to a ground station as displayed in the Figure 1(c). However, the UAVs cannot maintain a fixed position during their flight due to the minimum velocity constraint and the constraint of energy. As a result, they continuously revolve around the ground stations as closely as possible with a minimum velocity in order to preserve good communication channels.

According to the Figure 1(b), we note that the transmission power of the different UAVs is almost always set at its maximum level $P_{max}$. This is explained by the fact that, in the resulting solution of our algorithm, the UAVs are quite distant from each other during their communications with the ground stations. As a result, our resulting solution could be assimilated to the combination of two distant subsystems in which every one of them is composed of one UAV. However, The impact of the power transmission would be more tangible if we adopted a system which is composed of a large number UAVs. In that case, the communication channel would be highly impacted by the interference at the level of users as there are a lot of UAVs that are close to each other emitting the signals at the same time.

Remark: We can always consider working with a larger number of UAVs and ground stations. For instance, Figure 2 illustrates the optimal trajectories of UAVs in the case where we consider a system composed of 4 UAVs and 9 ground stations.

V-A Convergence of our proposed Algorithm

Let us study the convergence behavior of our algorithm. It can be observed from the Table I below that the max-min throughput increases in every iteration. We also note that our algorithm converges quickly to a final solution within about 10 iterations.

V-B Performance Analysis

Let us compare the performance of our proposed solution to a benchmark scheme with static access point, placed in the geometric center of the ground nodes. In fact, the baseline scheme allows us to achieve a performance of $R=17$ Bps/Hz, which is much more negligible to the performance achieved by our solution which is $R=208$ Bps/Hz. This result strongly highlight the benefit brought by UAV mobility in such problems.

Additionally, let us compare the performance of our proposed solution to the initial trajectory solution. We note that our proposed solution allows us to achieve results that are significantly better than those obtained from the initial solution. Indeed, the max min rate has increased from $R=70$ Bps/Hz to $R=208$ Bps/Hz for a duration $T=100$ s. This represents an improvement of nearly $150$ Bps/Hz in the overall performance.

We also note from the Table II that the UAVs have planned their movements so that the users have the same duration of communication with them. Additionally, we note that the ground stations have received roughly the same rate which is equal to $R=208$ Bps/Hz. Thus, we consider that the objective of achieving a fair performance among all users is achieved.

V-C Impact of the energy on the solution

In this section, we study the impact of the energy constraint on the performance of our system. So in the beginning, we denote $E_{m}$ as the energy consumed by the UAVs in the first example provided in Section V for $T=100s$. Next, we solve our problem for different values of $E_{max}$ and compare the results that we obtain with the first solution provided in Section V.

• •

$E_{max}$=0.9 $E_{m}$

We note that that the UAVs managed to find trajectories that are different from those we have in Figure 3. But, it succeeded to preserve the same max min rate $R=208$ bps/Hz.

• •

$E_{max}$= 0.6 $E_{m}$

Our algorithm succeeded to converge to a final solution which is represented in Figure 4. However, we note that the value of the max min rate has fallen drastically compared with the initial solution, from $R=208$ bps/Hz to $R=183.6$ bps/Hz.

• •

$E_{max}$= 0.3 $E_{m}$

We note that our algorithm failed to converge because it was not able to find a feasible set for our problem that is compliant with the different constraints of energy and velocity.

VI Conclusion

In this letter, we have investigated the max-min fairness problem between users in wireless communication with multiple UAVs. The main objective behind this problem was to design an efficient solution that jointly optimizes the UAVs’ trajectories, the transmission power and the use scheduling that meet certain quality of service requirements. To this end, we proposed an effective iterative algorithm that regroups the BCD and the SCA techniques. Then, we demonstrated the effectiveness and performance of our algorithm compared to the baseline circular path planning scheme.

Appendix A

In this section, we will solve (P3) using the successive convex optimization techniques. Let us first start by the constraints (3) and (8). The resulting sets of these constraints are not convex with respect to $\mathcal{V}$. Thereby, we start by introducing a slack variable $\lambda_{m}(n)$ and replacing it with $\|\mathbf{v}_{m}(n)\|$ . As a result, by defining $\Lambda=\{\lambda_{m}(n)$   ,$\forall m,n$} we obtain the following problem (P3.a) :

[TABLE]

$\Lambda$ is now considered as a new decision variable set. Additionally, we have also two new constraints (18) and (18). It is worth noting that at the optimal solution of this new problem (P3.a), we must have $\lambda_{m}(n)=\|\mathbf{v}_{m}(n)\|$. Otherwise, we can enlarge the feasible region corresponding to the constraint (18) by increasing $\lambda_{m}(n)$. Therefore, we can say with certainty that the problem (P3.a) is equivalent to the previous problem (P3).

In addition, with this new reformulation, the resulting set of the energy constraint (18) is now convex with respect to ($\mathcal{V,A},\Lambda$). The constraint (18) is also convex with respect to $\Lambda$. However, the new constraint (18) is non-convex with respect to $\mathcal{V}$.

Next, As stated earlier, we are trying to solve (P3.a) approximately by solving iteratively a series of convex problems that represent a lower bound to this latter. So, as we want to maximize $\mu$, our approach consists into shrinking the feasible set of the non-convex constraints by making some approximations in order to make them convex.

Specifically, we apply the successive convex optimization. So, we define the given velocity of an UAV in the $r^{th}$ iteration as $\mathcal{V}^{r}$ ={$\mathbf{v}^{r}_{m}(n)$   ,$\forall m,n\}$ which represents the solution of $\mathcal{V}^{r}$ computed at the $(r-1)^{th}$ iteration. Then, we use the following lemma:

Lemma 1**.**

If $f$ is a convex (respectively concave) function with respect to a certain variable $x\in X$. Then $f$ is lower (respectively upper) bounded by its first order Taylor expansion at any point $x,y\in X$. Thus:*

• •

$f(x)\geq f(y)+f’(y)(x-y)$,     if $f$ is convex

• •

$f(x)\leq f(y)+f’(y)(x-y)$,    if $f$ is concave

**

Since $f(x)$= $\|x\|^{2}$ is convex, we apply the previous lemma by taking $x=\mathbf{v}_{m}(n)$ and

$y=\mathbf{v}_{m}^{r}(n)$

Consequently, we have the following result:

[TABLE]

As a result, we define the new constraint

[TABLE]

This new constraint is convex with respect to ($\mathcal{V},\Lambda$).

Let us consider now the constraint (15) which is non-convex with respect to $\mathcal{B}$. In order to deal with it, we first write the left-hand side of the constraint as a difference between two concave functions with respect to the auxiliary variable $\mathcal{B}$. So, we have:

[TABLE]

where $\overline{R}_{k,m}(n)=\log_{2}\left(\sum_{\begin{subarray}{c}j=1\\ j\neq m\end{subarray}}^{M}B_{k,j}(n)+\sigma_{0}^{2}\right)$

As the first and the second terms of $R_{k,m}(n)$ are concave, we cannot study the convexity of $R_{k,m}(n)$. To tackle this problem, we can use the successive convex optimization in order to approximate $\bar{R}_{k,m}$(n) by a more controllable function at a given local point. Specifically, we define $\mathcal{B}^{r}$ ={$B^{r}_{k,m}(n)$   $\forall k,m,n$} which represents the given auxiliary variable calculated in the $(r-1)^{th}$ iteration.

Since $f(x)=\log(x)$ is concave, we apply the previous lemma by taking

$x=\sum_{\begin{subarray}{c}j=1\\ j\neq m\end{subarray}}^{M}B_{k,j}(n)+\sigma_{0}^{2}$ and $y=\sum_{\begin{subarray}{c}j=1\\ j\neq m\end{subarray}}^{M}B^{r}_{k,j}(n)+\sigma_{0}^{2}$

We obtain the following result:

$\overline{R}_{k,m}(n)\leq\overline{R}_{k,m}^{upper}(n)=\sum_{\begin{subarray}{c}j=1\\ j\neq m\end{subarray}}^{M}A_{k,j}(n)(B_{k,j}(n)-B_{k,j}^{r}(n))+C_{k,j}(n){}$

where the constants $A_{k,j}(n)$ and $C_{k,j}(n)$ are given by

[TABLE]

As a result, for a given trajectory $B^{r}$ defined in the $r^{th}$ iteration, the non-convex constraint (B-B) can be approximated by:

[TABLE]

Now, Let us consider the constraint (15). The resulting set of this constraint is not convex with respect to $\mathcal{Q}$. So, as explained in the beginning of this section, we need to find a lower bound to the right-hand side of this constraint which is convex.

As shown in (Appendix D), the function $B_{k,m,max}(n)=\frac{\beta_{0}P_{max}}{H^{2}+\|\mathbf{q}_{m}(n)-\mathbf{w}_{k}\|^{2}}{}$ is a concave surrogate function. As a result, after we define $\mathcal{Q}^{r}$ =$\{Q^{r}_{m}(n)$   ,$\forall m,n\}$ which represents the given trajectory calculated in the $(r-1)^{th}$ iteration, we have the following results:

[TABLE]

where the constants $D_{k,m}(n)$ and $F_{k,m}(n)$ are given by

[TABLE]

As a result, for a given trajectory $\mathcal{Q}^{r}$ defined in the $r^{th}$ iteration, the non-convex constraint (15) can be approximated by:

[TABLE]

Finally, for a given trajectory, velocity and the auxiliary variable ($\mathcal{Q}^{r},\mathcal{V}^{r},\mathcal{B}^{r}$), we obtain the problem (P4).

Thanks to the lower bounds that we adopted in the constraints (B-B), (15) and (18), the resulting set of all the constraints of the problem (P4) are convex. As a result, the problem (P4) is now considered as a convex optimization problem that can be solved efficiently using some predefined optimization solvers.

Appendix B Proof of the convergence and complexity of our Algorithm

B-A Proof of convergence

In order to prove that our algorithm converges, we first start by defining the objective values of our main problem (P2) $\mu$ as well as the sub-problems (P2.1) and (P4) at the $r^{th}$ iteration as $\mu^{r}(\mathcal{S})=\mu^{r}$ and $\mu^{lb,r}(\mathcal{A,Q,V,S,B})=\mu^{lb,r}$ , respectively.

Since (P2.1) is defined as a sub-problem of (P2) for fixed $(\mathcal{Q}^{r},\mathcal{B}^{r},\mathcal{V}^{r})$, then when we optimize (P2.1) over the user scheduling the objective value $\mu$ of the global problem (P2) is also maximized for a given $(\mathcal{Q}^{r},\mathcal{B}^{r},\mathcal{V}^{r})$. Then:

[TABLE]

Second, for given $(\mathcal{S}^{r+1},\mathcal{Q}^{r},\mathcal{V}^{r},\mathcal{B}^{r})$ in step 2.b) of our algorithm , we can express the following relationships:

[TABLE]

where the first equation (a) holds since the surrogate functions in (15) as well as the first-order Taylor expansions in (B-B) and (18) are tight at given local points. In other words, for given ($\mathcal{Q}^{r},\mathcal{V}^{r},\mathcal{B}^{r})$, we can consider that problem (P4) has the same objective value as that of problem (P2.2).

The second inequality (b) is derived from the fact that ($\mathcal{Q}^{r+1},\mathcal{V}^{r+1},\mathcal{B}^{r+1}$) represent the optimal solutions of problem (P4) for given ($\mathcal{S}^{r+1},\mathcal{Q}^{r},\mathcal{V}^{r},\mathcal{B}^{r}$) and thanks to the non-decreasing propriety of the objective function in this problem (because we are maximizing $\mu$), this inequality holds.

The third inequality (c) holds since the objective function of the approximation problem (P4) is considered as lower bound of that of the original sub-problem (P2.2).

Based on (19) and (20), we have the following result:

$\mu(\mathcal{S}^{r},\mathcal{Q}^{r},\mathcal{V}^{r},\mathcal{B}^{r})\leq\mu(\mathcal{S}^{r+1},\mathcal{Q}^{r+1},\mathcal{V}^{r+1},\mathcal{B}^{r+1})$

As a result, we can affirm that, in each iteration, the objective value $\mu$ of our main problem (P2) increases or at least remains unchanged. Additionally, since the objective value is upper bounded, the convergence of our algorithm is thus proved.

B-B Proof of complexity of our algorithm

In order to determine the complexity of our algorithm, we should determine the complexity of both subproblems (P2.1) and (P4) defined in section IV of our paper.

First, problem (P2.1) is a standard linear problem that depends only on the user-UAV association parameters. As explained in [12], the complexity of such problems is of order $\Theta\left(a^{2}b\right)$, where $a$ is the number of variables and $b$ is the number of constraints. Following this, it is easy to show that the complexity of (P2.1) is of order $\Theta\left((KMN+MN+KN+K)(KM)^{2}\right)$.

Second, the problem (P4) is defined as follows:

[TABLE]

Problem (P4) is a convex problem. As stated in [1], the complexity of (P4) is of order $\Theta\left(max(a^{3},a^{2}b,F)\right)$ where $a$ is the number of variables, $b$ is the number of constraints and $F$ is the cost of evaluating the first and second derivatives of the objective and constraint functions.

• •

First, it can be easily proved that $a=KMN+7MN+1$.

• •

Second, we can compute $b$ by summing the number of constraints in the equations (B-B)-(LABEL:a10) which are represented in the following table:

Thus, we can totally calculate the number of constraints as $b=KMN+4MN+2M(N-1)+K+M+\frac{M(M-1)N}{2}$.

• •

Since $F$ is the cost of evaluating the first and second derivatives of the objective and constraint functions, a simple approach consists into computing the first and second derivatives of the equations (B-B)-(LABEL:a10) with respect to each of the optimization variables ($\mu^{lb},\mathcal{B,Q,V,A},\Lambda$) separately and then summing the number of operations required to compute every single one of these quantities. By doing so, we get: $F\approx 122MN+(4M+9)KMN$.

It is clear that $a^{2}b$ is greater than $a^{3}$ and $F$. Thus, the complexity of (P4) is

[TABLE]

Finally, the complexity of one iteration of our algorithm is nothing but the sum of the complexities of the two subproblems which is

[TABLE]

If we assume that $K$, $M$ and $N$ are much larger than one ($\gg 1$), then the complexity of 1 iteration of our algorithm can be reduced to the order of

[TABLE]

Appendix C K-means clustering algorithm

C-A Definition

$K$-means clustering is an unsupervised learning algorithm which is mainly used in data mining and statistics. This iterative method aims to partition the data into $K$ clusters by allocating each point to the cluster with the nearest mean. In other words, this method clusters automatically similar data examples together.

The intuition behind $K$-means is an iterative procedure that is explained in the following algorithm:

  **Algorithm : K-means Clustering Algorithm **

Centroids initialization. 2. 2.

Repeat

1. (a)

Calculate distance between the points and the centroids. 2. (b)

Assign the points to the closest cluster. 3. (c)

Update the position of the centroids. 3. 3.

Until Convergence

As explained in the above algorithm, we first start by initializing the position of the centroids (points that represent the center of the clusters). In the second step, for each point, we calculate its distance to all the centroids and then assign it to the closest one. More formally, if we define $\mathcal{C}$ as the set of centroids, then each point $\mathbf{x}$ is attributed to a cluster based on:

$\arg\min\limits_{\mathbf{c}_{i}\in\mathcal{C}}(\|\mathbf{c}_{i}-\mathbf{x}\|)$

In the third step, we update the position of the centroids by taking the mean of the locations of all the points assigned to that centroid’s cluster. More precisely, the new positions of the centroids are given by :

[TABLE]

where $\mathcal{S}_{i}$ and $|\mathcal{S}_{i}|$ are respectively the set and the number of the ground stations that belongs to the $i^{th}$ cluster.

Finally, we repeat iteratively the second and third steps until there is no change in the position of centroids.

C-B Application of the $K$-means Clustering Algorithm

We apply the $K$-means clustering algorithm in order to divide the ground stations into $M$ clusters. First, we randomly initialize the position of the $M$ centroids. Then, as explained in the previous subsection, we first assign every ground station to the closest centroid and then we update the position of these centroids by taking the mean of the cluster. So, this process is repeated iteratively until convergence to a final solution.

Next, we attribute every cluster to a UAV $m$. So, every UAV $m$ will communicate only with the ground stations that are located in its given cluster.

The figures 5 and 6 illustrate the iterative process that we used. So, in the end, we obtain two clusters of ground stations in which the red dots are assigned to one UAV and the blue dots are assigned to the other UAV.

Afterwards, we provide a circular trajectory as an initialization scheme for every UAV in which the centroids that we obtained from the $K$-mean clustering algorithm will represent the center of the circular trajectory of every UAV. Furthermore, the radius of this trajectory circles will represent the average distance from the centroid to its assigned ground stations which belong to the same cluster. In other words the radius $r_{m}$ is defined by:

[TABLE]

where $\mathcal{S}_{m}$ and $|\mathcal{S}_{m}|$ represent respectively the set and the number of ground stations that belong to the cluster m.

Appendix D

Let us define $f(\mathbf{x})=\frac{1}{a\|\mathbf{x}\|^{2}+b}$ ,where $a$ and $b$ are positive constants. Then, the function $g$ defined by:

[TABLE]

is a concave surrogate function of $f$.

In order to prove this, we use the following lemma:

Lemma 2**.**

**

In order to prove that a function $g(\mathbf{x}|\mathbf{x}_{0})$ is a concave surrogate function of $f(\mathbf{x})$, $g(\mathbf{x}|\mathbf{x}_{0})$ needs to satisfy the following conditions

[TABLE]

**

Proof:

1) It is easy to say that $f(\mathbf{x})=g(\mathbf{x}|\mathbf{x}),\hskip 28.45274pt\forall\mathbf{x}$

2) We first calculate the gradient of $f$ at any given point. So, we have:

$\nabla f(\mathbf{x})=\frac{-2a}{(a\|\mathbf{x}\|^{2}+b)^{2}}\mathbf{x},\hskip 28.45274pt\forall\mathbf{x}$

Then we compute the gradient of $g$ with respect to $\mathbf{x}$. We obtain:

[TABLE]

Thus, it can be easily shown that:

          $\nabla g(\mathbf{x}|\mathbf{x})=\nabla f(\mathbf{x}),\hskip 2.84544pt\forall\mathbf{x}$

3) Let us show that: $g(\mathbf{x}|\mathbf{x}_{0})\leq f(\mathbf{x}),\hskip 2.84544pt\forall\mathbf{x}$

First we consider the mean value form of the Taylor series of the function $f(\mathbf{x})-g(\mathbf{x}|\mathbf{x}_{0})$ at given point $\mathbf{x}_{0}$. Then we have:

[TABLE]

where $\mathbf{\epsilon}$ is between $\mathbf{x}$ and $\mathbf{x}_{0}$.

Since we proved in 1) and 2) that $f(\mathbf{x})=g(\mathbf{x}|\mathbf{x})$ and $\nabla g(\mathbf{x}|\mathbf{x})=\nabla f(\mathbf{x})$, we can say that:

$f(\mathbf{x})-g(\mathbf{x}|\mathbf{x}_{0})=\frac{1}{2}(\mathbf{x}-\mathbf{x}_{0})^{T}(\nabla^{2}f(\mathbf{\epsilon})-\nabla^{2}g(\mathbf{\epsilon}|\mathbf{x}_{0}))(\mathbf{x}-\mathbf{x}_{0})$

Additionally, it can be easily shown that :

• •

$\nabla^{2}f(\mathbf{x})=\frac{-2a}{(a\|\mathbf{x}\|^{2}+b)^{2}}I+\frac{8a^{2}}{(a\|\mathbf{x}\|^{2}+b)^{3}}$

• •

$\nabla^{2}g(\mathbf{x}|\mathbf{x}_{0})=\frac{-2a}{b^{2}}I$

We note that the eigenvalues of $\nabla^{2}f(\mathbf{\epsilon})$ are greater than the eigenvalues of $\nabla^{2}g(\mathbf{\epsilon}|\mathbf{x}_{0})$. Thus $\nabla^{2}f(\mathbf{\epsilon})-\nabla^{2}g(\mathbf{\epsilon}|\mathbf{x}_{0})$ is a positive semi definite matrix.

As a result, $\frac{1}{2}(\mathbf{x}-\mathbf{x}_{0})^{T}(\nabla^{2}f(\mathbf{\epsilon})-\nabla^{2}g(\mathbf{\epsilon}|\mathbf{x}_{0}))(\mathbf{x}-\mathbf{x}_{0})\geq 0$

Consequently , $g(\mathbf{x}|\mathbf{x}_{0})\leq f(\mathbf{x}),\hskip 28.45274pt\forall\mathbf{x}$ ∎

Based on 1), 2) and 3), $g(\mathbf{x}|\mathbf{x}_{0})$ is a concave surrogate function of $f(\mathbf{x})$.

By substituting $\mathbf{x}=\mathbf{q}_{m}(n)-\mathbf{w}_{k}$, $\mathbf{x}_{0}=\mathbf{q}_{m}^{r}(n)-\mathbf{w}_{k}$, $b=H^{2}$, and $a=1$ and multiplying $f(\mathbf{x})$ and $g(\mathbf{x}|\mathbf{x}_{0})$ by $P_{max}\beta_{0}$ in the third condition of the lemma, we obtain the following result: $P_{max}\beta_{0}\left(-\frac{\|\mathbf{q}_{m}(n)-\mathbf{w}_{k}\|^{2}}{H^{4}}+D_{k,m}(n)(\mathbf{q}_{m}(n)-\mathbf{w}_{k})^{T}(\mathbf{q}_{m}^{r}(n)-\mathbf{w}_{k})+F_{k,m}(n)\right)\leq\frac{P_{max}\beta_{0}}{H^{2}+\|\mathbf{q}_{m}(t)-\mathbf{w}_{k}\|^{2}}$

where the constants $D_{k,m}(n)$ and $F_{k,m}(n)$ are given by:

[TABLE]