An Energy-Efficient Logistic Drone Routing Method considering Dynamic Drone Speed and Payload

1. Introduction

Over the past few decades, the development of technology, especially carbon fiber and lithium polymer batteries, has facilitated the use of drones for last mile delivery. For example, Zipline has completed over 14,000 medical supply deliveries in Rwanda since 2016 ([1]). In 2019, Google’s Wing and Amazon received Federal Aviation Administration (FAA) approval to begin commercial delivery by drone ([2,3]).

In recent years, researchers have shown a growing interest in last-mile delivery that incorporates drones ([4,5]). The drone delivery problems can be categorized as drone-only or truck-drone tandem problems. The first uses only drones for parcel delivery, while the second uses both drones and trucks to deliver. The alternative objectives of drone delivery problems are to minimize delivery time or minimize operational cost ([6]). Numerous variants of these problems have been suggested, such as heterogeneous drones ([7]) or multiple trucks ([8,9]).

Compared with traditional truck transportation, drone delivery is faster and less affected by road systems, in addition, it can save a lot of labor and time costs ([10]). These benefits allow logistics companies to dispatch drones for parcel delivery. On the other hand, a drone’s endurance is constrained by its battery’s capacity, which is further impacted by parcel weight, speed, and weather conditions ([11,12,13]). In the literature involving drone delivery, one typical assumption is that drones fly at constant or maximum speed, regardless of their payload or speed. Thus, speed is not taken into consideration during the decision-making procedure. In fact, operating drones to fly at a fixed speed may lead to increased energy consumption, thus resulting in an energy-infeasible route or an inability to serve long-distance customers because the energy consumption of drones is severely impacted by parcel weight and speed. Therefore, it is crucial to consider the impact of payload and speed on endurance in drone routing problems.

The motivation for considering flight speeds as decision variables is the influence of drone energy consumption on flight endurance and range. Figure 1(a) shows the power rate of drones under different payloads, using the model suggested by [12]. The circular marks indicate optimal drone speeds that maximize the drone endurance under any onboard battery capacity, which can be seen more clearly from Figure 1(c). Figure 1(b) shows the drone energy consumption per unit traveling distance with varying payloads. The triangular symbols indicate optimal drone speeds that maximize the total traveling distance with fixed battery capacity, which can be observed more distinctly from Figure 1(d). It can be seen from Figure 1(b) that the drone energy consumption per unit distance initially decreases and subsequently rises again as the speed increases. This is because there is a balance between the flight duration and the energy consumption. Flying at lower speeds requires less power, which leads to a longer flight duration. In contrast, higher speeds result in more power being consumed, but the makespan is reduced. Therefore, the flight range notably depends on its speed, and faster flight speeds are not always optimal. Consequently, higher speeds can cause faster deliveries but with reduced range and increased energy consumption. This is essential within the scope of drone delivery systems because, in some cases, the customer is sufficiently remote from the warehouse, and flying at maximum speed can lead to a shorter range that falls short of the distance necessary to reach the customer. It may be beneficial to fly at lower speeds which can improve range and serve more customers. Moreover, speed adjustments can enhance endurance, which enables drones to serve more customers. As such, it is imperative to factor in drone speeds during deliveries when considering energy consumption.

This paper introduces the multi-trip drone routing problem with variable flight speeds (DRP-VFS), an extension of the DRP proposed by [14]. In this problem, a fleet of homogeneous drones delivers multiple parcels to customers while minimizing energy consumption. In this case, drones can perform multiple trips and only be dispatched from and returned to the depot once per trip. While existing literature on the DRP assumes that drones fly at constant speeds in advance, the DRP-VFS considers drone speeds as decision variables when drones make deliveries.

This paper has the following contributions.

• We introduce and formally define a multi-trip drone routing problem, aiming to achieve operational cost minimization, where drone speeds are considered as decision variables rather than constants.
• We clearly consider drone energy consumption as a nonlinear function of both flight speed and payload rather than supposing that flight endurance and range are constants.
• We proposed a three-phased method integrating variable drone speed to solve DRP-VFS, which can improve the solution in terms of computational time to a greater extent than MILP.

The rest of this paper is structured as follows. In Section 2, we present a review of the literature. Section 3 presents the formal description of DRP-VFS and its MILP formulation. In Section 4, an efficient heuristic method proposed to solve large-scale instances is introduced. Section 5 demonstrates the case study and verifies the effectiveness of the proposed method. In the end, conclusions are drawn, and future research work is discussed in Section 6.

Figure 1. The required power, energy consumption per unit distance, endurance, and range of a drone versus speed and payload referenced from [15].

Figure 1. The required power, energy consumption per unit distance, endurance, and range of a drone versus speed and payload referenced from [15].

2. Literature Review

Over the past years, there have been increasing studies on the civil application of drones ([4,16]). These studies reported that operational planning for drone delivery is a well-studied problem. There already exists a lot of literature on the vehicle routing problem (VRP) involving conventional delivery vehicles and its variants ([17,18]). However, adapting the VRP to drones is faced with additional challenges since delivery by drones must account for extra constraints such as lower battery capacity and higher sensitivity of energy cost to payload and flight speed ([6]).

Most drone studies have assumed either a fixed endurance or a constant speed ([10,19,20,21,22,23,24]). Obviously, under this assumption, drones flying at higher speed have longer ranges. Other studies consider a constant range of UAVs. These two assumptions are equivalent when drone speeds are constant. However, both of them fail to consider the influence of payload and speed on the drone’s energy consumption, thus limiting the flexibility of drone delivery. For example, drones can extend their flight endurance by adjusting speed to serve customers further away from the depot.

Some studies have improved the fixed endurance/range model by considering drone power consumption dependent on payload while remaining unrelated to drone speed. [25] apply the UAV power rate function proposed by [26] and assume that the UAV maintains the maximum power during flight to calculate the load-dependent UAV speed during delivery. [27] propose a package delivery schedule utilizing drones and consider the drone’s power rate varies linearly with its payload. [6] propose two drone routing problems where customers are served by drones only. The authors consider the power rate function linearly dependent on payload weight and develop a simulated annealing (SA) algorithm to address the problems. The power rate function adopted in [6] pertains solely to hovering without encompassing that of actual flight. [14] adopt the same energy consumption as in [6], but instead of approximating the power rate function, they develop an exact algorithm adding logical and subgradient cuts to solve a DRP with time windows. [28] use the UAV power rate proposed by [29] as a piecewise linear function of package weight in their DRP. [30] also consider the influence of package weight on energy consumption using the same model as [6]. However, different from the proposed model based on varying speeds, the studies in the above papers only consider constant drone speeds. Moreover, the proposed energy-consumption models are mainly about hovering and do not consider the energy cost of actual flight status. Other studies apply an energy consumption model related to the payload and the drone speed during forward flight. [31] adopt an energy model by [26], which is linearly dependent on payload and speed. [32] also consider the energy consumption function introduced in [11] in the multiple flying sidekicks traveling salesman problem (mFSTSP). However, the UAV speeds are fixed in these two papers, limiting the potential time and energy savings from flying at varying speeds.

Some recent studies have considered variable flight speeds in drone routing problems. For instance, [33] considers the reciprocal of the flight speed to be approximately linear with the weight of loads. However, they consider drone speed as a function of the load, which can be regarded as a fixed parameter related to given loads rather than as a decision variable to be optimized. [34] consider drone speeds as decision variables and use the energy model introduced by [12]. Besides, [35] employs the multi-rotor UAV energy consumption model developed by [11] and considers power rate as a function of a UAV’s payload and speed. [35] investigate the mFSTSP with variable drone speeds. These two studies treat drone speeds as decision variables and demonstrate that operating drones at varying speeds contribute to delivery cost and time reductions, respectively. Moreover, [34] emphasizes that using fixed speeds may result in infeasibility in some instances. These two papers consider drone and truck tandem delivery, whereas our focus is solely on drone-only deliveries. As far as we know, this paper presents for the first time a DRP-VDS that considers UAV energy consumption and variable flight speeds.

Table 1 demonstrates a summary of relevant studies on the drone delivery problem reviewed in this paper. The studies are compared using the following factors: (i) the type of delivery (truck-drone tandem or drone-only), (ii) whether the endurance of drones is fixed or evaluated by energy consumption model, (iii) the drone speed (constant or variable), (iv) the objective function used, and (v) whether the solution methods are exact or heuristic. In the literature mentioned above, it is evident that only a few papers clearly consider the energy cost model and even fewer take drone speeds into account as decision variables. Because the finite drone range is susceptible to speed and payload weight, it is essential to consider variable drone speed and payload in drone delivery systems. In this study, a DRP-VFS is proposed to consider the energy constraints and variable speed and payload of the drone delivery systems for last-mile delivery systems. Also, a three-phased heuristic method is proposed to solve large-scale scenarios. The performance of the proposed heuristic is evaluated by comparing it to the exact method.

3. Model Formulation

This section presents a formulation of the DRP-VFS, and introduces piecewise linearization to approximate the nonlinear energy consumption function.

3.1. Problem Definition

The multi-trip drone routing problem with variable flight speeds (DRP-VFS) is defined on a complete directed graph $G=(N,A)$, where $N=\left\{0,1,\dots ,n+1\right\}$ represents the set of all nodes with the node 0 and $n+1$ corresponding to the depot (distribution center) from which all drones must originate and return. The set of customers is denoted by $N^{\prime }=\left\{1,\dots ,n\right\}$ associated with the distinct parcels. To simplify notation, we introduce $N^{+}=\left\{0,\dots ,n\right\}$ to indicate the set of nodes from which a drone can depart and $N^{-}=\left\{1,\dots ,n+1\right\}$ to describe the set of nodes that a drone may visit. Let $A=\left\{\left(i,j\right):i\in \left\{0\right\},j\in N^{\prime }\mathrm{and}i\in N^{\prime },j\in N^{-},i\ne j\right\}$ denote the set of arcs. Additionally, we define sets ${\delta }^{-}\left(i\right)$ and ${\delta }^{+}\left(i\right)$ to denote the forward and backward nodes of node i, respectively.

Each customer has a positive demand $q_i$ and a hard time window $\left[a_i,b_i\right]$. The drones are allowed to arrive early but need to wait until $a_i$ to begin serving the customer. A fleet of U homogeneous multi-rotor drones departs from the depot with a fully charged battery. Q represents the drone’s maximum payload, and we suppose that $q_i\le Q,\forall i\in N^{\prime }$, while E is the maximum energy provided by a fully charged battery. The service time for each customer is $s_i$. Every drone can implement multiple trips and serve several customers within one delivery. $d_{ij}$ represents flight distance on arc $\left(i,j\right)$, which is assumed to be symmetric. Since the drone speeds are decision variables in this paper, the flight time on arc $\left(i,j\right)$ is non-constant, as discussed later. The maximum speed limit is represented by $V_{max}$.

We make the following assumptions:

• Parcels are delivered by drones only.
• There is only one depot from which drones can depart.
• We neglect the time of loading parcels and swapping batteries.
• Without loss of generality, we ignore the influence of weather, i.e., wind impact is not considered.
• Drones can fly at a constant speed between two locations. The speed of each flight can vary.

3.2. Mathematical Model

DRP-VFS involves two sets of binary variables: $x_{ij}=1$ if a drone flies across arc $\left(i,j\right)$, and $x_{ij}=0$ otherwise. $z_{ij}=1$ if a trip ends in node i and is immediately succeeded by another trip with node j as the first customer. Additionally, there are six sets of continuous variables: $m_{ij}$ is the parcel’s weight transported across the arc $\left(i,j\right)$ (kg). $v_{ij}$ denotes the drone speed from node $i\in N^{+}$ to node $j\in N^{-}$ (m/s). ${\tau }_i$ is the service start time of customer $i\in N^{-}$ (second). $t_{ij}$ represents the time required to fly through arc $\left(i,j\right)$. $f_i$ denotes the total energy a drone consumes when it arrives at customer i (kWh). $e_{ij}$ indicates the energy consumed through arc $\left(i,j\right)$, measured in kWh.

The constraints of route feasibility can be expressed as follows:

$$\begin{array}{cc}\hfill \sum _{j\in {\delta }^{+}\left(i\right)}{x}_{ij}& =1\forall i\in N^{\prime },\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \sum _{j\in {\delta }^{-}\left(i\right)}{x}_{ji}& =1\forall i\in N^{\prime },\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \sum _{j\in N^{\prime }}{x}_{0j}& =\sum _{j\in N^{\prime }}{x}_{j,n+1}.\hfill \end{array}$$

(3)

Constraints (1) and (2) guarantee that each customer is served only once. Constraint (3) guarantees that the number of trips leaving the depot equals the number returning to the depot.

We apply demand constraints through

$$\begin{array}{cc}\hfill & \sum _{j\in {\delta }^{-}\left(i\right)}{m}_{ji}-\sum _{j\in {\delta }^{+}\left(i\right)}{m}_{ij}=q_i\forall i\in N^{\prime },\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill & m_{i,n+1}=0\forall i\in N^{\prime },\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill & m_{ij}\le Q{x}_{ij}\forall \left(i,j\right)\in A.\hfill \end{array}$$

(6)

Constraints (4) ensure that the demand of each customer can be satisfied, and also eliminate sub-tours. Consequently, the payload decreases to zero when the drone returns to the depot as specified by constraints (5). With constraints (6), the total weight of parcels carried by a drone in a trip must not exceed its payload capacity.

The following equations are formulated to enforce time-related constraints.

$$\begin{array}{cc}\hfill & {\tau }_i+s_i+d_{ij}/v_{ij}\le M\left(1-x_{ij}\right)+{\tau }_j\forall i\in N^{\prime },j\in N^{-},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & a_i\le {\tau }_i\le b_i\forall i\in N^{-},\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & {\tau }_i+s_i+\left(t_{i,n+1}+t_{0j}\right)\le M\left(1-z_{ij}\right)+{\tau }_j\forall i,j\in N^{\prime },i\ne j.\hfill \end{array}$$

(9)

Constraints (7) require that the total of the arrival time and service time in i and travel time across the arc $\left(i,j\right)$ not greater than the arrival time at the customer j (the immediate successor of i). Constraints (8) indicate that the arrival time to customer i falls within its designated time window $\left[a_i,b_i\right]$. Constraints (9) indicate the temporal relationship between successive trips executed by one drone, including the return time. Note that constraints (7) are nonlinear due to $\frac{1}{v_{ij}}$, we will introduce a piecewise linear method to convert them into linear ones in Section 3.3.

In this paper, we allow drones to change batteries and place parcels for new trips when they return to the depot. We adopt a 2-index formulation inspired by [14] to formulate the multi-trip characteristics of drone delivery. Therefore, we can achieve the reusability of drones by incorporating the following constraints:

$$\begin{array}{cc}\hfill \sum _{\begin{array}{c}i\in N^{\prime }\\ i\ne j\end{array}}{z}_{ij}& \le x_{0j}j\in N^{\prime },\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \sum _{\begin{array}{c}j\in N^{\prime }\\ j\ne i\end{array}}{z}_{ij}& \le x_{i,n+1}i\in N^{\prime },\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \sum _{j\in N^{\prime }}{x}_{0j}& -\sum _{i\in N^{\prime }}\sum _{\begin{array}{c}j\in N^{\prime }\\ j\ne i\end{array}}{z}_{ij}\le U.\hfill \end{array}$$

(12)

Constraints (10) and (11) establish a relationship between variables $\mathbf{x}$ and $\mathbf{z}$. Constraint (12) limits the number of drones that can be used.

To evaluate the energy usage during delivery, we adopt the power rate model according to [12].

$$P\left(v,w\right)=\frac{\sigma }{8}\rho sD{\mathsf{\Omega }}^2{R}^3(1+\frac{3v^2}{U_{tip}^2})+\frac{\kappa w^{\frac{3}{2}}}{\sqrt{2\rho D}}{\left(\sqrt{\frac{v^4}{4v_0^4}+1}-\frac{v^2}{2v_0^2}\right)}^{\frac{1}{2}}+\frac{1}{2}\rho S_{FP}{v}^3.$$

(13)

where $\sigma$ is the average profile drag coefficient, $\rho$ is the air density, s is the rotor solidity, $\mathsf{\Omega }$ is the angular velocity of the rotor, $U_{tip}$ is the blade tip speed, R is the rotor radius, D is the rotor disc area, $\kappa$ is the correction factor to ideal induced power, $v_0$ is the induced velocity in hovering, $S_{FP}$ is the equivalent flat plate area of the drone. w is the total weight of the drone, including the unloaded weight $m_v$ of the drone and the weight $m_p$ of parcels. Thus, $w=\left(m_v+m_p\right)g$ on the forward flight to customers, where g is the gravitational acceleration. To simplify the above formula, define $P_0=\frac{\sigma }{8}\rho sD{\mathsf{\Omega }}^2{R}^3$ and $P_i=\frac{\kappa w^{\frac{3}{2}}}{\sqrt{2\rho D}}$ as the blade profile power and induced power when hovering.

By transforming (13), we derive an energy consumption formula denoted as $E_0\left(v,m\right)$ ( $\mathrm{J}$/$\mathrm{m}$) for a drone’s travel over per unit distance at a constant speed v, and calculated as

$$E_0\left(v,w\right)=\frac{P\left(v,w\right)}{v}=P_0\left(\frac{1}{v}+\frac{3v}{U_{tip}^2}\right)+P_i{\left(\sqrt{v^{-4}+\frac{1}{4v_0^4}}-\frac{1}{2v_0^2}\right)}^{\frac{1}{2}}+\frac{1}{2}\rho S_{FP}{v}^2.$$

(14)

Thus, the energy consumption for drone on arc $\left(i,j\right)$, denoted as $e_{ij}$ can be expressed as

$$e_{ij}=\left\{\begin{array}{cc}{E}_0\left(v_{ij},w_{ij}\right)·d_{ij}\hfill & x_{ij}=1,\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(15)

where $d_{ij}$ and $w_{ij}=\left(m_v+m_{ij}\right)g$ are the distance and payload on arc $\left(i,j\right)$. It is worth noting that (15) is nonlinear when $x_{ij}=1$ and a piecewise linear method will be applied to deal with them in the next subsection.

Hence, drones’ energy constraints are written as

$$\begin{array}{cc}\hfill & f_0=0,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & f_i+e_{ij}\le M\left(1-x_{ij}\right)+f_j\forall \left(i,j\right)\in A,\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill & f_{n+1}\le E.\hfill \end{array}$$

(18)

Constraint (16) signifies that the initial energy consumption of each trip is set to zero. Constraints (17) indicate the energy consumed across the arc $\left(i,j\right)$. Constraint (18) guarantees the battery’s energy capacity constraint must be imposed.

The domain of the variables is shown below.

$$\begin{array}{cc}\hfill & x_{ij}\in \left\{0,1\right\}\forall \left(i,j\right)\in A,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & z_{ij}\in \left\{0,1\right\}\forall i,j\in N^{\prime },\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & m_{ij},e_{ij},t_{ij}\ge 0\forall \left(i,j\right)\in A,\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill & 0\le v_{ij}\le V_{max}\forall \left(i,j\right)\in A,\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & f_i\ge 0\forall i\in N,\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & {\tau }_i\ge 0\forall i\in N^{-}.\hfill \end{array}$$

(24)

3.3. Piecewise Linearization

Piecewise linearization can approximate a nonlinear univariate function using a set of non-negative variables of special ordered sets type 2 (SOS2), where at most two adjacent variables can be positive, while the sum of all variables equals 1 ([36]). For univariate functions like $1/v$, the piecewise linearization approximation is realized by sampling K points $\left[V_1,V_2,\dots ,V_K\right]$ as breakpoints on the v-axis. After that, the linear segments $\left[\left(V_1,1/V_1\right),\dots ,\left(V_K,1/V_K\right)\right]$ are used to approximate $1/v$.

In this case, we define a set of continuous variable ${\alpha }_{ij}^k$ for each breakpoint. The value of $v_{ij}$ on arc $\left(i,j\right)$ is represented by

$$v_{ij}=\sum _{k=1}^K{V}_k{\alpha }_{ij}^k.$$

(25)

The approximation value of the reciprocal of speed $1/v_{ij}$ is $1/{\tilde{v}}_{ij}$ which can be calculated by (26).

$$\frac{1}{{\tilde{v}}_{ij}}=\sum _{k=1}^K\frac{{\alpha }_{ij}^k}{V_k}.$$

(26)

The SOS2 variables ${\alpha }_{ij}^k$ need to satisfy the following constraints.

$$\begin{array}{cc}\hfill & \sum _{k=1}^K{\alpha }_{ij}^k=1,\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill & 0\le {\alpha }_{ij}^k\le 1k=1,2,\dots ,K.\hfill \end{array}$$

(28)

Most modern MIP solvers, such as Gurobi and CPLEX, can address specially ordered sets of types 1 and 2 automatically. Consequently, constraints (25)–(27) are sufficient to produce the correct computation. Consequently, the nonlinear constraints in (7) are approximated by replacing $1/v_{ij}$ with $1/{\tilde{v}}_{ij}$.

Similarly, we introduce a set of SOS3 variables ${\beta }_{ij}^{st}$ to linearly approximate the two-variable energy consumption function dependent on drone speed and parcel weight in (14). Similar to the SOS2, the SOS3 variables require that at most three adjacent components be greater than 0. The difference is that the set of SOS2 variables is one-dimensional, while the set of SOS3 variables is two-dimensional. Consider again S sampling coordinates $\left[V_1,V_2,\dots ,V_S\right]$ on the v-axis and T sampling coordinates $\left[M_1,M_2,\dots ,M_T\right]$ on the m-axis. The function $E_0\left(v,m\right)$ is determined by each breakpoint $\left(V_s,M_t\right)\left(s=1,\dots ,S;t=1,\dots ,T\right)$.

In our problem, the value of $v_{ij}$ and $m_{ij}$ on arc $\left(i,j\right)$ are represented by

$$\begin{array}{cc}\hfill & v_{ij}=\sum _{s=1}^S{V}_s{\beta }_{ij}^{st},\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill & m_{ij}=\sum _{t=1}^T{M}_t{\beta }_{ij}^{st}.\hfill \end{array}$$

(30)

The approximation value of the energy consumption is ${\tilde{e}}_{ij}$ which can be calculated by using (31).

$${\tilde{e}}_{ij}=\sum _{s=1}^S\sum _{t=1}^T{E}_0\left(V_s,M_t\right){\beta }_{ij}^{st}.$$

(31)

The SOS3 variables ${\beta }_{ij}^{st}$ should satisfy (32) and (33).

$$\begin{array}{cc}& \sum _{s=1}^S\sum _{t=1}^T{\beta }_{ij}^{st}=1,\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill & 0\le {\beta }_{ij}^{st}\le 1s=1,2,\dots ,S,t=1,2,\dots ,T.\hfill \end{array}$$

(33)

However, unlike SOS1 and SOS2, modern MIP solvers don’t have an automatic syntax for imposing SOS3. For this reason, we utilize the triangle method proposed by [37] to impose corresponding constraints. Consider the rectangle corresponding to intervals $\left[V_s,V_{s+1}\right)$ and $\left[M_t,M_{t+1}\right)$ and the two triangles generated by its diagonal $\left[\left(V_s,M_t\right)\left(V_{s+1},M_{t+1}\right)\right]$. We introduce two sets of binary variables $u_{ij}^{st}$ and $l_{ij}^{st}$ to represent the upper and lower triangles in the rectangle, respectively, with dummy values $u_{ij}^{0*}=u_{ij}^{*0}=u_{ij}^{S*}=u_{ij}^{*T}=0$ and $l_{ij}^{0*}=l_{ij}^{*0}=l_{ij}^{S*}=l_{ij}^{*T}=0$ at the boundaries. The additional constraints are shown as follow:

$$\begin{array}{cc}\hfill & \sum _{s=1}^{S-1}\sum _{t=1}^{T-1}\left(u_{ij}^{st}+l_{ij}^{st}\right)=1,\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill & {\beta }_{ij}^{st}\le u_{ij}^{st}+l_{ij}^{st}+u_{ij}^{s,t-1}+l_{ij}^{s-1,t-1}+u_{ij}^{s-1,t-1}+l_{ij}^{s-1,t}.\hfill \end{array}$$

(35)

The constraint (34) requires that just one triangle be chosen from all possible triangles. Constraints (35) guarantee that the only ${\beta }_{ij}^{st}$ values greater than 0 can be those corresponding to the three vertices of such a triangle. Accordingly, the nonlinear energy-related constraints in (17) are approximated by replacing $e_{ij}$ with ${\tilde{e}}_{ij}$.

The objective function is minimizing the total energy consumption during delivery.

$$\begin{array}{cc}\hfill min& \sum _{\left(i,j\right)\in A}{e}_{ij},\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \left(1\right)-\left(35\right).\hfill \end{array}$$

(36)

The model described above is a mixed-integer linear programming (MILP) model which can be solved efficiently by commercial MIP solvers such as Gurobi. While MILP solvers find optimal solutions, the computational time required is substantial, even for scenarios with several customers. In Section 4, we propose a three-phased heuristic algorithm aimed at finding suboptimal solutions within computational time constraints.

4. Solution Method

The DRP-VFS extends the classical VRP and is an NP-hard problem. Owing to the NP-hardness of the DRP-VFS, a three-phased heuristic is proposed to effectively tackle practical large-scale scenarios.

The solution of DRP-VFS can be represented by $\mathcal{R}=\left\{{\mathbf{r}}_1,{\mathbf{r}}_2,\dots ,{\mathbf{r}}_{\mathcal{K}}\right\}$, where $\mathcal{K}$ is the number of trips. The ${\mathbf{r}}_k$ contains a vector of customers that a drone will visit, which can be expressed as ${\mathbf{r}}_k=\left\{r_1,r_2,\dots ,r_{C_k}\right\}$, where $C_k$ is the number of customers assigned to the kth trip.

4.1. Initialization

There are numerous heuristic approaches that can efficiently search for a feasible solution for the VRP, such as CW savings algorithm and sweep algorithm. In our method, we first apply an improved k-means algorithm that includes capacity and Euclidean distance as constraints to partition customers into several subsets (Algorithm ). Next, a traveling salesman problem with time windows (TSPTW) is solved within each subset to determine the customer sequence a drone serves. This algorithm is executed while satisfying the capacity, energy, and time window constraints. The number of clusters is calculated by the total demands and capacity of the drone as

$$\mathcal{K}=\lceil \sum _{i=1}^n{q}_i/Q\rceil$$

(37)

The $\mathcal{K}$ centroids are randomly initialized in the customers’ coordinate range. The Euclidean distances between each customer and all $\mathcal{K}$ centroids are calculated. Sort customers by demand from largest to smallest and assign each customer to the nearest centroids. The total demands for each subset of clusters cannot exceed the drone’s capacity. The selected customer will be assigned to the next nearest centroid if the capacity constraint is not satisfied. The centroid of kth cluster $\left(x_k,y_k\right)$ is calculated based on its members using (38).

$$\begin{array}{c}\hfill x_k=\sum _{i=1}^{C_k}{x}_i/C_k,\\ \hfill y_k=\sum _{i=1}^{C_k}{y}_i/C_k.\end{array}$$

(38)

The iterative procedure is repeated until the cluster is unchanged or the maximum iteration is reached.

The TSPTW is a special case of DRP-VFS since only one drone is used. Thus, we just need to remove the reusability constraints (10)-(12) and modify constraints (3) to the following.

$$\sum _{r\in {\mathbf{r}}_k}{x}_{0r}=\sum _{r\in {\mathbf{r}}_k}{x}_{r,n+1}=1k=1,2,\dots ,\mathcal{K}.$$

(39)

Therefore, an initial feasible solution can be obtained quickly.

4.2. Local Search

In this section, we propose a speed optimization model to adjust the flight speed in a trip and calculate the cost of a DRP-VFS solution, and the Adaptive Large Neighborhood Search (ALNS) algorithm is used to enhance the solution quality.

4.2.1. Speed Optimization

The speed optimization problem (SOP) involves adjusting the flight speed in order to minimize energy consumption. The SOP is defined on a route $\mathbf{r}=\left\{r_0,r_1,\dots ,r_{C+1}\right\}$ consisting of customers served by a single drone, where $r_0$ and $r_{C+1}$ represent the depot corresponding to 0 and $n+1$, respectively. $v_i$ and $m_i$ denote the flight speed and payload carried by a drone between node $r_i$ and its successor. The SOP is formulated as follows:

$$min\sum _{i=0}^C{e}_i,$$

(40)

subject to

$$\begin{array}{cc}\hfill & m_i=\sum _{j=i+1}^C{q}_{j}i=0,\dots ,C,\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill & m_0\le Q,\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill & {\tau }_i+s_i+d_i/v_i\le {\tau }_{i+1}i=0,\dots ,C,\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill & a_i\le {\tau }_i\le b_{i}i=1,\dots ,C,\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill & e_i=E_0\left(v_i,m_i\right)d_{i}i=0,\dots ,C,\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill & f_0=0,\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill & f_i+e_i\le f_{i+1}i=0,\dots ,C,\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill & f_{C+1}\le E.\hfill \end{array}$$

(48)

The nonlinear energy function is linearized using SOS2 variables introduced in Section 3.3. Since $m_i$ is calculated in advance using (42), the $E_0$ between node $r_i$ and its successor depends only on $v_i$. We introduce ${\gamma }_i^s$ to linearly approximate the energy consumption function. Hence, the value of $v_i$ and the approximate value of ${\tilde{e}}_i$ can be represented by

$$\begin{array}{cc}\hfill & v_i=\sum _{s=1}^S{V}_s{\gamma }_i^s,\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill & {\tilde{e}}_i=\sum _{s=1}^S{E}_0\left(V_s,m_i\right){\gamma }_i^s.\hfill \end{array}$$

(50)

To deal with $1/v_i$, the same method can be used. However, to speed up the solving process, we introduce a convexification method, i.e., replace $1/v_i$ with an auxiliary variable ${\lambda }_i$ to avoid reintroducing integer variables:

$${\tau }_i+s_i+d_i{\lambda }_i\le {\tau }_{i+1}.$$

(51)

Meanwhile, the following inequality should be imposed to ensure ${\lambda }_i$ converges to $1/v_i$.

$${\lambda }_i\ge 1/v_i.$$

(52)

It can be seen that constraint (52) is a quadratic and convex constraint. Therefore, the SOP is modeled as a mixed-integer quadratically constrained programming (MIQCP) and can be solved by commercial solvers.

4.2.2. ALNS heuristics

In this subsection, an adaptive large neighborhood search algorithm is adopted to improve the initial solution obtained by Section 4.1. The ALNS framework has demonstrated successful results in solving multiple variants of vehicle routing problems ([38]). By removing selected nodes from the solution and reinserting them, we obtain a neighborhood of the given solution. The operators are dynamically chosen based on their previous performance. For this reason, each operator is given a score that increases if it enhances the current solution. A new solution is accepted if it meets the criteria specified by the local search scheme, such as simulated annealing.

In this paper, we apply two removal operators and two insertion operators to effectively destroy and repair a DRP-VFS solution. In addition, for every $\Gamma$ iteration, the TSPTW model is solved once within each trip to adjust the customer sequence. The framework of the ALNS algorithm is presented in Algorithm . Figure 2(a)-(c) depict the initial solution, the solution after applying two destroy operators and the solution after executing two repair operators.

The destroy operators involve primarily removing several nodes from the current solution and placing them in a removal list $\mathcal{L}$. We present the removal operators involved in our approach as follows:

(1)

Random removal: This operator removes several nodes at random from the current solution.

(2)

Worst removal: This operator removes the highest-cost node from the solution, where the cost is determined by solving SOP.

The insertion operators are applied to repair a partly destroyed solution by reinserting the nodes from $\mathcal{L}$ back into the solution. The insertion operators employed in our algorithm include:

(1)

Greedy insertion: This operator repeatedly removes a node from $\mathcal{L}$ and inserts it into the lowest cost position of a route.

(2)

Regret insertion: An obvious disadvantage of greedy insertion is that it defers node insertion to later iterations where few feasible moves are available. The regret operator in our algorithm uses a 2-regret criteria. Define $\Delta c_i^j$ as the cost change incurred by inserting node i into the route where the cost is jth cheapest. The 2-regret criteria inserts the node i based on $i^{*}=\underset{i\in \mathcal{L}}{\mathrm{argmax}}\left(\Delta c_i^2-\Delta c_i^1\right)$, where $\Delta c_i^1$ and $\Delta c_i^2$ are the best and second-best insertion of node i. We iterate the procedure until no more nodes in $\mathcal{L}$ can be inserted.

We describe the adaptive weight adjustment process below. The operators are chosen by the roulette-wheel method. In the beginning, all operators are equally likely. They are updated as $P_i^{t+1}=P_i^t\left(1-\mathsf{\Theta }\right)+\mathsf{\Theta }{\pi }_i/{\omega }_i$ during program execution, where $P_i^t$ is the probability of operator i at the tth iteration, $\mathsf{\Theta }$ is the parameter of roulette wheel, ${\pi }_i$ denotes the score of operator i and ${\omega }_i$ indicates times it was selected in the previous $\Gamma$ iterations. The score reflects the performance of each operator. When a new best solution is found, the related operators’ scores are increased by ${\sigma }_1$. If the new solution outperforms the current one, the score is increased by ${\sigma }_2$. If the new solution falls short of the current one but meets the acceptance criteria, the score is increased by ${\sigma }_3$.

Our ALNS algorithm applies the simulated annealing as an acceptance criterion. We calculate the probability in line to determine whether to update the current solution and decrease the initial temperature T according to $T=\mathsf{\Phi }T$, where $\mathsf{\Phi }$ is a constant indicating the cooling rate. The algorithm produces the best solution that has been found after $\mathsf{\Lambda }$ iterations.

4.3. Assignment

In this section, we describe a list schedule algorithm to minimize the number of drones. The delivery time, denoted by h, is the latest allowable departure time from the depot, whereas the arrival time, denoted by l, indicates when the drone reaches the depot upon finishing the delivery. The two times are calculated by using (53).

$$\begin{array}{cc}\hfill h& ={\tau }_1-d_0/v_0,\hfill \\ \hfill l& ={\tau }_C+s_C+d_C/v_C,\hfill \end{array}$$

(53)

where ${\tau }_1$ and ${\tau }_C$ are the starts of service time at the first and last nodes in a route, $d_0,v_0$ and $d_C,v_C$ is the distance and speed at the first and last flights in a route, and $s_C$ is the service time of the last node.

We allocate $\mathcal{K}$ routes to U drones using the list scheduling algorithm (Algorithm ). First, sort routes by delivery time in ascending order. A binary heap is used to get the route with minimum arrival time quickly. Next, traverse the $\mathcal{K}$ routes; if the delivery time of ${\mathbf{r}}_k$ is greater or equal to the arrival time of the heap top route, then ${\mathbf{r}}_k$ is appended to this route. Otherwise, a new drone should be added. The time complexity is $\mathcal{O}\left(\mathcal{K}log\mathcal{K}\right)$. In conclusion, our heuristic method finds a near-optimal solution ${\mathcal{R}}^{*}$ for DRP-VFS.

5. Results and Discussion

In this section, numerical experiments have been carried out to verify the computational efficiency of the proposed heuristic method for the DRP-VFS. Firstly, the computational performance of the ALNS algorithm is revealed by comparing its solution quality and CPU time to those of the exact method with small-scale scenarios in Section 5.2. Secondly, the performance of ALNS implementation in large-scale instances is demonstrated in Section 5.3. Finally, we compare the solutions generated by fixed-speed and variable-speed in Section 5.4. It is demonstrated that the energy consumption of variable-speed solutions decreases notably compared to that of fixed-speed solutions. All computational work was conducted on the computer configuration with an Intel i7 CPU and 16 GB RAM. The heuristic algorithm was coded in Python 3.7 and Gurobi 9.5.0 was used to solve the mathematical models.

5.1. Parameter Settings

We use two sets of benchmark instances introduced by [14], named $SetA$ and $SetB$. The $SetA$ is created based on [39] and [6], which contains 10-50 customers. The $SetB$ is an extension of Solomon’s instance, which contains 10-100 customers. We use $SetA$ to validate the effectiveness of MILP and ALNS implementation and energy savings with variable flight speeds. We use $SetB$ to validate the efficiency of ALNS implementation on large-scale instances. We execute the MILP implementation once for each instance. We execute the ALNS algorithm 10 times for each instance and calculate the average objective value and computational time. The parameter settings of drones and the ALNS algorithm are shown in Table 2 and Table 3. For $SetA$ instances, the battery energy capacity is set to $E=0.27$ kWh; For $SetB$ instances, we set $E=0.027$ kWh.

5.2. Performance Comparison between MILP and ALNS Implementations

In this section, we compare the performance of the MILP implementation with the ALNS implementation to demonstrate that the ALNS algorithm is effective in finding near-optimal solutions to small DRP-VFS instances. The results of MILP and ALNS are presented in Table 4. CPU is the computational time taken to solve the problem. The last three columns in Table 4 are the percentage differences of ALNS versus MILP in terms of energy consumption, total flight distance, and computational time. The reason why ALNS needs more drones than MILP is that the time window is not fully utilized when calculating the delivery time and arrival time using (53). For example, when appending ${\mathbf{r}}_1$ to ${\mathbf{r}}_2$, the delivery time of ${\mathbf{r}}_2$ is possible to postpone backward while satisfying the time window constraint. Nevertheless, ALNS does not take this scenario into account due to the challenges associated with its implementation.

We can observe that the CPU time of MILP grows exponentially with the number of customers. For instance, the average CPU time increases from $144.2$ $\mathrm{s}$ to $666.9$ $\mathrm{s}$, up by $362.5$ when solving the 20-customer problem compared to the 15-customer problem. The average CPU time of the ALNS increases from $123.5$ $\mathrm{s}$ to $149.4$ $\mathrm{s}$, a growth of $21.3$ . Moreover, for the first two instances with 10 customers, ALNS finds the optimal solution. It can be obviously found from Table 4 that the difference in energy consumption and travel distance between ALNS and MILP is less than 10 , which reveals that ALNS is capable of finding near-optimal solutions in small scenarios.

5.3. Performance of ALNS on Large-Scale Instances

Additionally, the ALNS algorithm is used to solve large-scale instances. Table 5 indicates the results of $SetB$ instances with 100 customers. Trips is the number of trips that drones would execute. The ALNS algorithm can find feasible solutions for large instances within a limited time. However, due to its high spatial complexity, the MILP model can’t be solved. As shown in Table 5, the average computational time is $689.6$ $\mathrm{s}$, indicating the high computational efficiency of the proposed heuristic algorithm.

5.4. Comparison between Variable-Speed and Fixed-Speed

This section provides the energy saving achievable by enabling drones to fly at variable speeds. The MILP implementation is run with varying speeds and three fixed speeds (10/20/30 $\mathrm{m}$/$\mathrm{s}$) to generate solutions for the DRP-VFS. Table 6 shows that the average total energy consumption can be reduced by up to $46.61$ against the 10 $\mathrm{m}$/$\mathrm{s}$ case, $1.03$ against the 20 $\mathrm{m}$/$\mathrm{s}$ case and $15.13$ against the 30 $\mathrm{m}$/$\mathrm{s}$ case. The drone speeds in variable speed cases are all located at maximum-range speed since this is the most energy-efficient speed for drones. The results in Table 6 demonstrate that variable flight speeds result in lower energy consumption than fixed speeds in most instances. Reducing energy consumption could enable the utilization of smaller-capacity batteries, resulting in weight savings and possibly extending the flight range.

6. Conclusion and Future Work

This paper extends the multi-trip DRP with time windows to incorporate variable flight speeds. We proposed a MILP formulation for DRP-VFS that minimizes total cost while considering the energy consumption model, flight speed, payload weight, and drone reuse. To tackle practical instances, we proposed a three-phased ALNS heuristic algorithm. Numerical results indicate that optimizing drone speeds is an important consideration for drone delivery. Variable-speed resulted in about 46 and 15 improvements in total energy consumption compared to fixed-speed 10 $\mathrm{m}$/$\mathrm{s}$ and 30 $\mathrm{m}$/$\mathrm{s}$ cases, respectively. Numerical experiments also revealed that the three-phased heuristic performed well in both solution quality and CPU times in small and large DRP-VFS instances.

There are multiple opportunities for further research in this scope. For instance, exact solution methods can be explored to address the nonlinear energy function. To extend the flight range of drones, multiple depots or recharging stations might be investigated. More destruction and repair operations can be applied to improve solution quality. Additionally, the problem could be raised by considering heterogeneous drones or collaborative delivery of drones and trucks.