Drone-Fleet-Enabled Logistics: A Joint Design of Flight Trajectory and Package Delivery

Abstract

In this work, we focus on a drone-fleet-enabled package delivery scenario, in which multiple drones fly from a start point and cooperatively deliver packages to the ground users in the presence of a number of no-fly zones (NFZs). We first mathematically model the package delivery scenario in a rigorous manner. Then, a package value maximization problem is established to optimize the flight trajectory and package delivery under the constraints of drone load and collision as well as NFZs. The formulated problem is a highly challenging mixed-integer non-convex problem. To facilitate solving it, we transform the formulated problem into an equivalent problem with special structure by using some appropriate transformations, based on which a low-complexity algorithm with favorable performance is developed using the penalty convex–concave procedure method. Finally, numerical results demonstrate the superiority of the proposed solution.

1. Introduction

Over the past few years, drones have attracted extensive attention in the fields of wireless communications, the internet of things, logistics, and so on, due to their high mobility, autonomy, and flexibility [1]. In a wide variety of drone applications, the flight trajectory design of the drone is a pivotal but challenging task, since it is required not only to prevent collisions between drones, but also to avoid no-fly zones (NFZs), such as military bases, government agencies, and airports [2]. To tackle this challenging task, some representative research efforts have emerged in recent years, such as [3,4,5,6,7]. Specifically, the authors of [3,4] studied the problem of minimizing drone power consumption through the joint design of wireless resource management and path planning under NFZ constraints in drone-enabled communication systems. In [5], by carefully designing the drone trajectory in the presence of multiple NFZs, the authors maximized the uplink sum rate among all ground users. Note that only a single-drone scenario was examined in [3,4,5]. For a multi-drone scenario (i.e., drone fleet), the authors of [6,7] studied secure drone-enabled communication systems constrained by NFZs and collision avoidance, in which the system confidentiality was maximized by carefully planning the drone fleet trajectory.

We remark that the drones in the above literature [3,4,5,6,7] were only used as flying base stations to enhance the performance of the existing terrestrial wireless networks. However, as pointed out by [8], the use of drones for package delivery is increasingly becoming a new industry trend. In the drone delivery scenario, drones may not be able to deliver packages for all users at the same time due to load constraints, which therefore calls for a joint design of the flight trajectory and package delivery. Recently, only several works such as [9,10] have focused on the trajectory design in a drone delivery scenario, with the goal of minimizing the propulsion energy consumed by a single drone [9] or the flight path length of a drone fleet [10]. However, in [9,10], the authors implicitly assumed that there were no NFZs in the network, and did not take into account drone load limits, which obviously is not applicable for the real world. As a result, when considering the NFZ and load constraints, how to jointly design the flight trajectory and package delivery of drones has still not been tackled appropriately, which motivates this work.

Our main contributions can be summarized as follows. First, we conduct rigorous mathematical modeling of the drone-enabled package delivery scenario, in which a fleet of drones takes off from a start point and cooperatively delivers packages to users on the ground in the presence of multiple NFZs. Then, we establish a package value maximization problem by jointly optimizing the flight trajectory and package delivery under the constraints of drone load and collision as well as NFZs. The formulated problem is a highly challenging mixed-integer non-convex problem. To facilitate solving it, we transform the formulated problem into an equivalent problem with special structure by using some appropriate transformations. On this basis, we propose a low-complexity algorithm with favorable performance based on the penalty convex–concave procedure method. Finally, numerical results demonstrate the superiority of the proposed solution.

2. System Model

Let us consider a package delivery scenario using a drone fleet, as shown in Figure 1, in which there are U drones flying from the start point to the end point, aiming to deliver packages to K ground users within T time slots. For ease of illustration, let us define $\mathcal{U}\triangleq \{1,2,\cdots ,U\}$, $\mathcal{K}\triangleq \{1,2,\cdots ,K\}$, and $\mathcal{T}\triangleq \{1,2,\cdots ,T\}$ to be the drone set, the user set, and the time slot set, respectively. Let ${\mathbf{c}}_u\left[t\right]\in {\mathbb{R}}^{2\times 1}$ represent the horizontal location of the u-th drone in slot $t\in \mathcal{T}$. Then, to avoid collision between drones, we should have

$$\begin{array}{cc}\hfill & \parallel {\mathbf{c}}_u\left[t\right]-{\mathbf{c}}_{\tilde{u}}\left[t\right]\parallel \ge L^{min},t\in \mathcal{T},u\in \mathcal{U},\tilde{u}\in \mathcal{U},\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & {\mathbf{c}}_u\left[1\right]={\mathbf{c}}^{\mathrm{str}},{\mathbf{c}}_u\left[T\right]={\mathbf{c}}^{\mathrm{end}},\hfill \end{array}$$

(2)

where $L^{min}$ is the minimum distance between drones in each slot, and ${\mathbf{c}}^{\mathrm{str}}\in {\mathbb{R}}^{2\times 1}$ and ${\mathbf{c}}^{\mathrm{end}}\in {\mathbb{R}}^{2\times 1}$ are the horizontal locations of the start and end points, respectively. In addition, to bound the flight speed of drone u, we further have

$$\begin{array}{c}\hfill \parallel {\mathbf{c}}_u\left[t\right]-{\mathbf{c}}_u[t-1]\parallel \le L_u^{max},u\in \mathcal{U},t\in \mathcal{T}\backslash \left\{1\right\},\end{array}$$

(3)

where $L_u^{max}$ is the u-th drone’s maximum traveling distance in a slot. We define $\mathbf{C}\triangleq {\left({\mathbf{C}}_u\right)}_{u\in \mathcal{U}}$ with ${\mathbf{C}}_u\triangleq {\left({\mathbf{c}}_u\left[t\right]\right)}_{t\in \mathcal{T}}$ as the trajectory design of the drone fleet.

We consider that there are K packages in the depot at the start point, and each package belongs to a different user. For ease of illustration, we refer to the package belonging to user k as package k. We assume that package k can be divided into $I_k\ge 1$ sub-packages, denoted by set ${\mathcal{I}}_k\triangleq \{1,2,\cdots ,I_k\}$. Each sub-package $i\in {\mathcal{I}}_k$ has a weight of $w_{i,k}$ (in kilos) and thus the total weight of package k is given by $w_k\triangleq {\sum }_{i\in {\mathcal{I}}_k}{w}_{i,k}$. Note that each sub-package can be delivered by one drone at most and each package can be delivered cooperatively by multiple drones. Let $W_u$ (in kilos) denote the load limit of the u-th drone. Since different drones may have different load limits, we need to select a suitable number of drones for each package. Specifically, let $d_{i,k,u}\in \{0,1\}$ represent the delivery indicator variable of sub-package $i\in {\mathcal{I}}_k$, where $d_{i,k,u}=1$ denotes that sub-package $i\in {\mathcal{I}}_k$ is delivered by drone u, and $d_{i,k,u}=0$, otherwise. Then, we have the following constraints:

$$\begin{array}{cc}\hfill & d_{i,k,u}\in \{0,1\},i\in {\mathcal{I}}_k,k\in \mathcal{K},u\in \mathcal{U},\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill & \sum _{u\in \mathcal{U}}{d}_{i,k,u}\le 1,i\in {\mathcal{I}}_k,k\in \mathcal{K},\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill & \sum _{u\in \mathcal{U}}\sum _{i\in {\mathcal{I}}_k}{d}_{i,k,u}\in \{0,I_k\},k\in \mathcal{K},\hfill \end{array}$$

(6)

where the constraint in (5) is due to the fact that each sub-package can only be delivered by one drone at most and the constraint in (6) comes from our stipulation that package k is delivered only when all sub-packages of package k are delivered. In addition, to prevent overload of drone u, the following constraint should be met:

$$\begin{array}{cc}\hfill & \sum _{k\in \mathcal{K}}\sum _{i\in {\mathcal{I}}_k}{d}_{i,k,u}{w}_{i,k}\le W_u,u\in \mathcal{U}.\hfill \end{array}$$

(7)

We define $\mathbf{D}\triangleq {\left(d_{i,k,u}\right)}_{i\in {\mathcal{I}}_k,k\in \mathcal{K},u\in \mathcal{U}}$ as the package delivery design. Let $s_k\ge 0$ denote the value of the k-th package. Then, the total value of the packages delivered by the fleet can be expressed as $S\left(\mathbf{D}\right)\triangleq {\sum }_{k\in \mathcal{K}}\frac{s_k}{I_k}{\sum }_{u\in \mathcal{U}}{\sum }_{i\in {\mathcal{I}}_k}{d}_{i,k,u}$. Note that in terms of the physical meaning, the value $s_k$ can represent the importance of package k or the priority of user k.

In order to deliver package k, the drone, such as the u-th drone, must fly to the delivery zone around user k, as depicted in Figure 1. We determine the delivery zone of user k according to whether drone u can detect the uplink signal of user k. As such, it is required to calculate the received signal-to-noise ratio at drone u during flight. As in [11], we assume that the channels from the users to the drones are only dominated by line-of-sight links. Then, we can compute the channel power gain from user k to drone u at slot t as $g_{k,u}\left[t\right]\triangleq \frac{g_0}{h^2+{\parallel {\mathbf{c}}_u\left[t\right]-{\mathbf{a}}_k\parallel }^2}$, where $g_0$ denotes the channel power gain at a reference distance of one meter, ${\mathbf{a}}_k\in {\mathbb{R}}^{2\times 1}$ is the horizontal location of user k, and h is a constant, representing the altitude of each drone [1]. As a result, at drone u, the received signal-to-noise ratio from user k is given by ${\gamma }_{k,u}\left[t\right]=\frac{p_k{g}_{k,u}\left[t\right]}{N_0}$. Here, $p_k$ is the transmission power of user k and $N_0$ is the noise power. In view of the above notations, we now formally determine that the delivery zone of user k is the area where ${\gamma }_{k,u}\left[t\right]$ is greater than or equal to a pre-defined threshold ${\theta }_u$. Here, ${\theta }_u$ is related to the signal detection capability of drone u. Sub-package $i\in {\mathcal{I}}_k$ can be delivered only when drone u arrives at the delivery zone of user k, so we have

$$\begin{array}{cc}\hfill & d_{i,k,u}\le \sum _{t\in \mathcal{T}}𝟙[{\gamma }_{k,u}\left[t\right]\ge {\theta }_u],i\in {\mathcal{I}}_k,k\in \mathcal{K},u\in \mathcal{U},\hfill \end{array}$$

(8)

where $𝟙[·]$ represents an indicator function.

Finally, it should be noted that the drones are prohibited from crossing the NFZs when delivering packages, as illustrated in Figure 1. We assume that there exist N NFZs in the network, denoted by $\mathcal{N}\triangleq \{1,2,\cdots ,N\}$, and characterize the n-th NFZ as a cylindrical volume with a height greater than or equal to h [6]. We define ${\mathbf{b}}_n\in {\mathbb{R}}^{2\times 1}$ and $Z_n$ to be the horizontal coordinate center and radius of NFZ n, respectively. Then, to keep away from each NFZ, the following constraint should be satisfied:

$$\begin{array}{c}\hfill \parallel {\mathbf{c}}_u\left[t\right]-{\mathbf{b}}_n\parallel \ge Z_n,n\in \mathcal{N},u\in \mathcal{U},t\in \mathcal{T}.\end{array}$$

(9)

3. Problem Establishment

In this work, we would like to maximize the total value of the packages delivered by the fleet, i.e., $S\left(\mathbf{D}\right)$, through jointly optimizing the trajectory design $\mathbf{C}$ and the package delivery design $\mathbf{D}$ subject to the constraints in (1)–(9). Specifically, we have the following optimization problem.

$$\begin{array}{cc}\hfill \underset{\mathbf{C},\mathbf{D}}{max}& S\left(\mathbf{D}\right)\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \left(1\right)–\left(9\right).\hfill \end{array}$$

(10)

The problem in (10) is a challenging mixed-integer non-convex problem with the indicator function. To solve it in a computationally tractable manner, some necessary transformations are required, as detailed below.

To be specific, we first eliminate the indicator function in (8). Let $M_{k,u}\left[t\right]$ denote a large constant greater than any possible value of $\parallel {\mathbf{c}}_u\left[t\right]-{\mathbf{a}}_k{\parallel }^2$. Then, by utilizing the big-M approach and introducing a new binary variable $\mathbf{X}\triangleq {\left(x_{k,u}\left[t\right]\right)}_{k\in \mathcal{K},u\in \mathcal{U},t\in \mathcal{T}}$, we can equivalently convert the constraint in (8) to

$$\begin{array}{cc}\hfill & x_{k,u}\left[t\right]\in \{0,1\},k\in \mathcal{K},u\in \mathcal{U},t\in \mathcal{T},\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & d_{i,k,u}\le \sum _{t\in \mathcal{T}}{x}_{k,u}\left[t\right],i\in {\mathcal{I}}_k,k\in \mathcal{K},u\in \mathcal{U},\hfill \end{array}$$

(12)

$$\begin{array}{c}\parallel {\mathbf{c}}_u\left[t\right]-{\mathbf{a}}_k{\parallel }^2\le M_{k,u}\left[t\right](1-x_{k,u}\left[t\right])\hfill \\ +\frac{g_0{p}_k}{\theta N_0}-h^2,k\in \mathcal{K},u\in \mathcal{U},t\in \mathcal{T}.\hfill \end{array}$$

(13)

Second, by introducing a new binary variable $\mathbf{B}\triangleq {\left(b_k\right)}_{k\in \mathcal{K}}$, we equivalently convert the constraint in (6) to

$$\begin{array}{cc}\hfill & b_k\in \{0,1\},k\in \mathcal{K},\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & b_k=\frac{1}{I_k}\sum _{u\in \mathcal{U}}\sum _{i\in {\mathcal{I}}_k}{d}_{i,k,u},k\in \mathcal{K}.\hfill \end{array}$$

(15)

Third, we eliminate the binary variables $\mathbf{D}$, $\mathbf{X}$, and $\mathbf{B}$. Using the similar method adopted in [12], we can equivalently transform the binary constraints in (4), (11) and (14) to the following continuous constraints:

$$\begin{array}{cc}\hfill & d_{i,k,u}\in [0,1],i\in {\mathcal{I}}_k,k\in \mathcal{K},u\in \mathcal{U},\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & x_{k,u}\left[t\right]\in [0,1],k\in \mathcal{K},u\in \mathcal{U},t\in \mathcal{T},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill & b_k\in [0,1],k\in \mathcal{K},\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill & d_{i,k,u}(1-d_{i,k,u})\le 0,i\in {\mathcal{I}}_k,k\in \mathcal{K},u\in \mathcal{U},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & x_{k,u}\left[t\right](1-x_{k,u}\left[t\right])\le 0,k\in \mathcal{K},u\in \mathcal{U},t\in \mathcal{T},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & b_k(1-b_k)\le 0,k\in \mathcal{K}.\hfill \end{array}$$

(21)

Last, based on the above transformations, the problem in (10) can be equivalently transformed to the following problem.

$$\begin{array}{cc}\hfill & \underset{\mathbf{C},\mathbf{D},\mathbf{X},\mathbf{B}}{max}S\left(\mathbf{D}\right)\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \left(1\right)–\left(3\right),\left(5\right),\left(7\right),\left(9\right),\left(12\right),\left(13\right),\left(15\right)–\left(21\right).\hfill \end{array}$$

(22)

Note that, the constraints in (1), (9) and (19)–(21) are non-convex. Fortunately, the left-hand sides of these constraints can be written as the difference of two convex functions. In the following section, we will see that this special structure greatly facilitates solving the problem in (22).

4. Low-Complexity Solution

Due to the special structure mentioned above, the problem in (22) can be regarded as a difference-of-convex problem. As such, we propose a low-complexity algorithm for solving the problem in (22) based on the penalty convex–concave procedure, which generally consists of two steps, namely relaxation and linearization. The overall procedure for solving the problem in (22) is shown in Figure 2.

4.1. Relaxation

In this step, by adding slack variables to (1), (9) and (19)–(21), and penalizing the sum of the corresponding violations, we obtain a penalized difference-of-convex problem as follows.

$$\begin{array}{cc}\hfill & \underset{\mathbf{C},\mathbf{D},\mathbf{X},\mathbf{B},\epsilon ,\sigma ,\alpha ,\beta ,\zeta }{max}S\left(\mathbf{D}\right)-\xi P(\mathit{\epsilon },\mathit{\sigma },\mathit{\alpha },\mathit{\beta },\mathit{\zeta })\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \left(2\right),\left(3\right),\left(5\right),\left(7\right),\left(12\right),\left(13\right),\left(15\right)–\left(18\right),\hfill \\ \hfill & -\parallel {\mathbf{c}}_u\left[t\right]-{\mathbf{c}}_{\tilde{u}}{\left[t\right]\parallel }^2\le -{\left(L^{min}\right)}^2\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & +{\epsilon }_{u,\tilde{u}}\left[t\right],u\in \mathcal{U},\tilde{u}\in \mathcal{U},t\in \mathcal{T},\hfill \\ \hfill & -\parallel {\mathbf{c}}_u\left[t\right]-{\mathbf{b}}_n{\parallel }^2\le -Z_n^2\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill & +{\sigma }_{n,u}\left[t\right],n\in \mathcal{N},u\in \mathcal{U},t\in \mathcal{T},\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill & d_{i,k,u}(1-d_{i,k,u})\le {\alpha }_{i,k,u},i\in {\mathcal{I}}_k,k\in \mathcal{K},u\in \mathcal{U},\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & x_{k,u}\left[t\right](1-x_{k,u}\left[t\right])\le {\beta }_{k,u}\left[t\right],k\in \mathcal{K},u\in \mathcal{U},t\in \mathcal{T},\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill & b_k(1-b_k)\le {\zeta }_k,k\in \mathcal{K},\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill & {\epsilon }_{u,\tilde{u}}\left[t\right]\ge 0,u\in \mathcal{U},\tilde{u}\in \mathcal{U},t\in \mathcal{T},\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill & {\sigma }_{n,u}\left[t\right]\ge 0,n\in \mathcal{N},u\in \mathcal{U},t\in \mathcal{T},\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill & {\alpha }_{i,k,u}\ge 0,i\in {\mathcal{I}}_k,k\in \mathcal{K},u\in \mathcal{U},\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill & {\beta }_{k,u}\left[t\right]\ge 0,k\in \mathcal{K},u\in \mathcal{U},t\in \mathcal{T},\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill & {\zeta }_k\ge 0,k\in \mathcal{K}.\hfill \end{array}$$

(33)

Here,

$$\begin{array}{c}\hfill P(\mathit{\epsilon },\mathit{\sigma },\mathit{\alpha },\mathit{\beta },\mathit{\zeta })\triangleq \sum _{t\in \mathcal{T}}\sum _{u\in \mathcal{U}}\sum _{\tilde{u}\in \mathcal{U}}{\epsilon }_{u,\tilde{u}}\left[t\right]+\sum _{t\in \mathcal{T}}\sum _{n\in \mathcal{N}}\sum _{u\in \mathcal{U}}{\sigma }_{n,u}\left[t\right]+\\ \hfill \sum _{i\in {\mathcal{I}}_k}\sum _{k\in \mathcal{K}}\sum _{u\in \mathcal{U}}{\alpha }_{i,k,u}+\sum _{t\in \mathcal{T}}\sum _{k\in \mathcal{K}}\sum _{u\in \mathcal{U}}{\beta }_{k,u}\left[t\right]+\sum _{k\in \mathcal{K}}{\zeta }_k,\end{array}$$

denotes the total violation, where $\mathit{\epsilon }\triangleq {\left({\epsilon }_{u,\tilde{u}}\left[t\right]\right)}_{u\in \mathcal{U},\tilde{u}\in \mathcal{U},t\in \mathcal{T}}$, $\mathit{\sigma }\triangleq {\left({\sigma }_{n,u}\left[t\right]\right)}_{n\in \mathcal{N},u\in \mathcal{U},t\in \mathcal{T}}$, $\mathit{\alpha }\triangleq {\left({\alpha }_{i,k,u}\right)}_{i\in {\mathcal{I}}_k,k\in \mathcal{K},u\in \mathcal{U}}$, $\mathit{\beta }\triangleq {\left({\beta }_{k,u}\left[t\right]\right)}_{k\in \mathcal{K},u\in \mathcal{U},t\in \mathcal{T}}$, and $\mathit{\zeta }\triangleq {\left({\zeta }_k\right)}_{k\in \mathcal{K}}$ are the slack variables for the constraints in (1), (9) and (19)–(21), respectively. The coefficient $\xi >0$ in the problem in (23) represents a penalty parameter and the increase in $\xi$ will force the total violation $P(\mathit{\epsilon },\mathit{\sigma },\mathit{\alpha },\mathit{\beta },\mathit{\zeta })$ to zero. In particular, if $P(\mathit{\epsilon },\mathit{\sigma },\mathit{\alpha },\mathit{\beta },\mathit{\zeta })=0$, the problems in (22) and (23) are equivalent.

4.2. Linearization

In this step, we linearize the concave terms in (24)–(28); that is, $-\parallel {\mathbf{c}}_u\left[t\right]-{\mathbf{c}}_{\tilde{u}}{\left[t\right]\parallel }^2$, $-\parallel \mathbf{c}\left[t\right]-{\mathbf{b}}_n{\parallel }^2$, $-{\left(d_{i,k,u}\right)}^2$, $-{\left(x_{k,u}\left[t\right]\right)}^2$, and $-{\left(b_k\right)}^2$, by applying the first-order Taylor expansion at the given point $({\mathbf{C}}^{\left(j\right)},{\mathbf{D}}^{\left(j\right)},{\mathbf{X}}^{\left(j\right)},{\mathbf{B}}^{\left(j\right)})$, where ${\mathbf{C}}^{\left(j\right)}\triangleq {\left({\mathbf{c}}_u^{\left(j\right)}\left[t\right]\right)}_{u\in \mathcal{U},t\in \mathcal{T}}$, ${\mathbf{D}}^{\left(j\right)}\triangleq {\left(d_{i,k,u}^{\left(j\right)}\right)}_{i\in {\mathcal{I}}_k,k\in \mathcal{K},u\in \mathcal{U}}$, ${\mathbf{X}}^{\left(j\right)}\triangleq {\left(x_{k,u}^{\left(j\right)}\left[t\right]\right)}_{k\in \mathcal{K},u\in \mathcal{U},t\in \mathcal{T}}$, ${\mathbf{B}}^{\left(j\right)}\triangleq {\left(b_k^{\left(j\right)}\right)}_{k\in \mathcal{K}}$, and $j=1,2,\cdots$ denotes an iteration index. As such, at iteration j, the constraints in (24)–(28) can be approximated as the following linear constraints:

$$\begin{array}{c}\parallel {\mathbf{c}}_u^{\left(j\right)}\left[t\right]-{\mathbf{c}}_{\tilde{u}}^{\left(j\right)}{\left[t\right]\parallel }^2-2{({\mathbf{c}}_u^{\left(j\right)}\left[t\right]-{\mathbf{c}}_{\tilde{u}}^{\left(j\right)}\left[t\right])}^{\mathrm{T}}({\mathbf{c}}_u\left[t\right]-{\mathbf{c}}_{\tilde{u}}\left[t\right])\hfill \\ \le -{\left(L^{min}\right)}^2+{\epsilon }_{u,\tilde{u}}\left[t\right],u\in \mathcal{U},\tilde{u}\in \mathcal{U},t\in \mathcal{T},\hfill \end{array}$$

(34)

$$\begin{array}{c}-\parallel {\mathbf{c}}_u^{\left(j\right)}\left[t\right]-{\mathbf{b}}_n{\parallel }^2-2{({\mathbf{c}}_u^{\left(j\right)}\left[t\right]-{\mathbf{b}}_n)}^{\mathrm{T}}({\mathbf{c}}_u\left[t\right]-{\mathbf{c}}_u^{\left(j\right)}\left[t\right])\hfill \\ \le -Z_n^2+{\sigma }_{n,u}\left[t\right],n\in \mathcal{N},u\in \mathcal{U},t\in \mathcal{T},\hfill \end{array}$$

(35)

$$\begin{array}{c}{\left(d_{i,k,u}^{\left(j\right)}\right)}^2+d_{i,k,u}(1-2d_{i,k,u}^{\left(j\right)})\hfill \\ \le {\alpha }_{i,k,u},i\in {\mathcal{I}}_k,k\in \mathcal{K},u\in \mathcal{U},\hfill \end{array}$$

(36)

$$\begin{array}{c}{\left(x_{k,u}^{\left(j\right)}\left[t\right]\right)}^2+x_{k,u}\left[t\right](1-2x_{k,u}^{\left(j\right)}\left[t\right])\hfill \\ \le {\beta }_{k,u}\left[t\right],k\in \mathcal{K},u\in \mathcal{U},t\in \mathcal{T},\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill & {\left(b_k^{\left(j\right)}\right)}^2+b_k(1-2b_k^{\left(j\right)})\le {\zeta }_k,k\in \mathcal{K}.\hfill \end{array}$$

(38)

Here, ${(·)}^{\mathrm{T}}$ represents the transpose operator. By replacing (24)–(28) with (34)–(38), we can obtain a convex approximation subproblem at iteration j, as detailed below.

$$\begin{array}{cc}\hfill & \underset{\mathbf{C},\mathbf{D},\mathbf{X},\mathbf{B},\epsilon ,\sigma ,\alpha ,\beta ,\zeta }{max}S\left(\mathbf{D}\right)-{\xi }^{\left(j\right)}P(\mathit{\epsilon },\mathit{\sigma },\mathit{\alpha },\mathit{\beta },\mathit{\zeta })\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \left(2\right),\left(3\right),\left(5\right),\left(7\right),\left(12\right),\left(13\right),\left(29\right)–\left(38\right).\hfill \end{array}$$

(39)

Here, ${\xi }^{\left(j\right)}$ is the penalty coefficient at iteration j. The problem in (39) is a convex optimization problem, which, therefore, can be solved optimally using existing software toolkits such as CVX.

4.3. Algorithm Summary

The detailed algorithm for solving the problem in (22) is summarized in Algorithm 1. Algorithm 1 begins with an arbitrary initial point $({\mathbf{C}}^{\left(0\right)},{\mathbf{D}}^{\left(0\right)},{\mathbf{X}}^{\left(0\right)},{\mathbf{B}}^{\left(0\right)})$ and a small penalty coefficient ${\xi }^{\left(0\right)}$, and then increases the penalty coefficient with a constant $\nu >1$ at each iteration until a limit bound ${\xi }_{max}$ is achieved to ensure $P(\mathit{\epsilon },\mathit{\sigma },\mathit{\alpha },\mathit{\beta },\mathit{\zeta })\mathit{\to }\mathit{\to }0$. It is worth mentioning that Algorithm 1 is able to converge to a stationary point of the problem in (22). The proof is similar to that in [13], which is omitted here for the sake of brevity. The complexity of Algorithm 1 is dominated by solving the problem in (39). If this problem is solved via CVX, the complexity is given by $\mathcal{O}(A^{3.5}log(1/ϵ))$ [14], where $A=UT(N+U+2K+2)+2K(U+1)$ and $ϵ>0$ represent the number of optimization variables and the given solution accuracy, respectively.

Algorithm 1: Low-complexity Algorithm for the problem in (22)

1: Choose an arbitrary initial point $({\mathbf{C}}^{\left(0\right)},{\mathbf{D}}^{\left(0\right)},{\mathbf{X}}^{\left(0\right)},{\mathbf{B}}^{\left(0\right)})$, ${\xi }^{\left(0\right)}>0$, ${\xi }_{max}$, and $\nu >1$. 2: Set $j\leftarrow 0$. 3: repeat 4:    Solve the problem in (39) to obtain an optimal solution $({\mathbf{C}}^{(j+1)},{\mathbf{D}}^{(j+1)},{\mathbf{X}}^{(j+1)},{\mathbf{B}}^{(j+1)})$. 5:    Set the optimal value of the problem in (39) to ${\widehat{V}}^{(j+1)}$. 6:    Set ${\xi }^{(j+1)}\leftarrow min\{\nu {\xi }^{\left(j\right)},{\xi }_{max}\}$. 7:    Set $j\leftarrow j+1$. 8: until$\frac{{\widehat{V}}^{(j+1)}-{\widehat{V}}^{\left(j\right)}}{{\widehat{V}}^{\left(j\right)}}\le 1\times {10}^{-5}$ or $j>50$.

5. Numerical Results

To facilitate demonstrating the performance of our proposed scheme, we consider that the users and NFZs are randomly distributed in a 2D area of $1.2\times 1.2$ ${\mathrm{km}}^2$. Unless otherwise mentioned, the simulation settings are given as follows: ${\mathbf{c}}^{\mathrm{str}}={(-550,-550)}^{\mathrm{T}}$, ${\mathbf{c}}^{\mathrm{end}}={(550,550)}^{\mathrm{T}}$, $N=4$, $I_k=1$, $p_k=0.01$ Watt for all $k\in \mathcal{K}$, $W_u=55$ kilos, ${\theta }_u=-1.5$ dB, $L_u^{max}=50$ meters for all $u\in \mathcal{U}$, $h=100$, $L^{min}=15$ meters, $g_0={10}^{-6}$ dB, $N_0={10}^{-12}$ Watt, $Z_n=150$ meters if n is odd and $Z_n=100$ meters, otherwise. In addition, we consider that the weight and the value of package $k\in \mathcal{K}$ follow a uniform distribution in the range $[10,40]$ and each sub-package of package k has the same weight. For Algorithm 1, we randomly choose $\nu$ from $(1.0,10]$ in each iteration and set ${\xi }^{\left(0\right)}=1\times {10}^{-6}$ and ${\xi }_{max}=1\times {10}^6$.

Figure 3 shows the flight trajectory of the drone fleet, where we set $T=35$, $K=5$, $U=3$, and $L^{min}=50$ m. From Figure 3, some observations can be made as follows. First, we can observe that each drone in the fleet does not fly into the NFZs, demonstrating that the proposed Algorithm 1 is able to effectively avoid the NFZs. Second, from their trajectory curves, we can easily verify that there are no collisions between drones. (It is noted that although the trajectory curves of drones are overlapped at some points, they are not overlapped at the same slot. For instance, for the first cross-point of drones 1 and 2, the horizontal coordinates of the two drones are ${\mathbf{c}}_1\left[8\right]$ and ${\mathbf{c}}_2\left[7\right]$, and for the second cross-point, the horizontal coordinates are ${\mathbf{c}}_1\left[24\right]$ and ${\mathbf{c}}_2\left[25\right]$, where ${\mathbf{c}}_1\left[8\right]={\mathbf{c}}_2\left[7\right]$ and ${\mathbf{c}}_1\left[24\right]={\mathbf{c}}_2\left[25\right]$. We can see that drones 1 and 2 actually do not meet at the same slot and thus, they will not collide.) Last, we can find that some drones (e.g., drones 1 and 3) will touch the users’ delivery zone, since they are required to perform package delivery tasks.

Figure 4 and Figure 5 compare our proposed scheme with some representative baselines at $T=100$ and $K=10$, where the results are averaged over 500 trials. To be specific, under the maximum load constraint in (7) the following baselines are considered. For Baseline 1, multiple packages with the highest value are delivered by the fleet. For Baseline 2, multiple lightest packages are delivered by the fleet. For Baseline 3, multiple packages closest to the end point are delivered by the fleet. Note that for each baseline, the flight trajectory of each drone can be achieved via Algorithm 1.

More specifically, Figure 4 and Figure 5 plot the total value $S\left(\mathbf{D}\right)$ versus the number of drones U and the load limit $W_u$, respectively. From Figure 4 and Figure 5, it first can be seen that the performance of all schemes improves as U ($W_u$) increases, due to the increase in fleet carrying capacity. However, with the increase in U ($W_u$), the marginal performance gain of adding more drones (storage capacity) into the fleet decreases, because the overall carrying capacity of the fleet is gradually large enough to carry all users’ packages. In addition, we also see that Baseline 2 significantly outperforms (slightly underperforms) Baseline 1 and Baseline 3 when $W_u$ is relatively small (large), e.g., $W_u\ge 40$ ($W_u<40$), indicating that the value- and distance-based schemes cannot work well in the area of each drone having a relatively small load limit. Finally, we can observe that the proposed scheme is superior to all baselines, because it can adapt well to the changes in the key system parameters through carefully optimizing the package delivery and fleet trajectory.

6. Discussion

Compared to the existing studies, the major innovativeness of our work lies in the following two aspects. The first point is that we focus on a drone-fleet-enabled package delivery scenario, in which multiple drones fly from a start point and cooperatively deliver packages to the ground users in the presence of a number of NFZs. This scenario is rarely considered in the current literature, but it is rigorously modeled mathematically in our paper. The second point is that by formulating a mixed-integer non-convex optimization problem and utilizing some appropriate transformation techniques (such as eliminating the indicator function and binary variables), we successfully develop a low-complexity algorithm based on the penalty convex–concave procedure method and achieve a joint design of drone fleet trajectory and package delivery.

The advantages of the joint design of drone fleet trajectory and package delivery proposed in this work lie in the following two aspects. On the one hand, under our design, the drone fleet can successfully deliver all packages of the ground users neither passing through the existing NFZs nor colliding with each other. On the other hand, the performance of our proposed joint design is better than the representative baselines due to its good adaptability to different system design parameters, such as the number of drones and the drone load limits, etc. These advantages suggest that the proposed joint design in this paper can address the challenges presented in Section 1 (i.e., the Introduction), and we believe that it can provide a theoretical reference for drone-enabled package delivery in real-world settings.

There are several interesting research directions worthy of extending the results of this paper in the future. (i) Energy-efficient drone-fleet-enabled package delivery: In this paper, we aim to maximize the total value of the packages delivered under the constraints of NFZs and the drone load, as well as the collision avoidance. We have not considered the energy consumption of drones. However, due to the limited capacity of the onboard battery of drones, the energy consumption of the drones needs to be considered when designing the flight trajectories of drones [15,16,17]. Thus, it is necessary to maximize the total value of delivered packages under the constraint of the energy consumption of the drone fleet. (ii) Optimal 3D trajectory of drone fleet: In this paper, we ignore the changes in the drone height in order to simplify the flight trajectory design of the drone fleet. However, each drone cannot only change its horizontal coordinate, but also its height; that is to say, the flying environment of the drone is 3D [18,19,20]. Therefore, in our future work, we are going to consider more realistic scenarios to optimize the 3D flight trajectory and package delivery design to maximize the total value of delivered packages.

7. Conclusions

In this paper, we studied a package delivery scenario supported by a drone fleet, in which a drone fleet flies from the start point in the presence of multiple NFZs and delivers packages to the ground users in a cooperative manner. We first rigorously modeled the package delivery scenario. Then, a package value maximization problem was established to optimize the flight trajectory and package delivery under the constraints of drone load and collision as well as NFZs. The formulated problem was a challenging mixed-integer non-convex problem. To facilitate solving it, we transformed the formulated problem into an equivalent problem with special structure by using some appropriate transformations. On this basis, we devised a low-complexity algorithm with promising performance by using the penalty convex–concave procedure method. Numerical results finally demonstrated the superiority of the proposed solution.