# Anomalous thermal transport across the superionic transition in ice

Rong Qiu  Qiyu Zeng  Department of Physics, National University of Defense
Technology, Changsha 410073, China  Hunan Key Laboratory of Extreme Matter and
Applications, National University of Defense Technology, Changsha 410073,
China  Han Wang  Laboratory of Computational Physics, Institute of Applied
Physics and Computational Mathematics, Beijing 100088, China  Dongdong Kang
Department of Physics, National University of Defense Technology, Changsha
410073, China  Hunan Key Laboratory of Extreme Matter and Applications,
National University of Defense Technology, Changsha 410073, China  Xiaoxiang
Yu<EMAIL_ADDRESS>Jiayu Dai<EMAIL_ADDRESS>Department of Physics,
National University of Defense Technology, Changsha 410073, China  Hunan Key
Laboratory of Extreme Matter and Applications, National University of Defense
Technology, Changsha 410073, China

###### Abstract

Superionic ices with highly mobile protons within the stable oxygen sub-
lattice occupy an important proportion of the phase diagram of ice and widely
exist in the interior of icy giants and throughout the universe. Understanding
the thermal transport in superionic ice is vital for the thermal evolution of
icy planets. However, it is highly challenging due to the extreme
thermodynamic conditions and dynamical nature of protons, beyond the
capability of the traditional lattice dynamics and empirical potential
molecular dynamics approaches. In this work, by utilizing the deep potential
molecular dynamics approach, we investigate the thermal conductivity of ice-
VII and superionic ice-VII” along the isobar of $p=30\ \rm{GPa}$. A non-
monotonic trend of thermal conductivity with elevated temperature is observed.
Through heat flux decomposition and trajectory-based spectra analysis, we show
that the thermally-activated proton diffusion in ice-VII and superionic ice-
VII” contribute significantly to heat convection, while the broadening in
vibrational energy peaks and significant softening of transverse acoustic
branches lead to a reduction in heat conduction. The competition between
proton diffusion and phonon scattering results in anomalous thermal transport
across the superionic transition in ice. This work unravels the important role
of proton diffusion in the thermal transport of high-pressure ice. Our
approach provides new insights into modeling the thermal transport and
atomistic dynamics in superionic materials.

As one of the most abundant substances in Earth and the universe, ice is of
vital importance from a scientific perspective and attracts wide research
interests. Especially, in the interior conditions of icy moons, where pressure
ranges from 2 GPa to hundreds of GPa and temperature ranges from 300 K to 4000
K, high-pressure ice phases (VII/VII”/X) is expected to widely exist [1].
These phases present the same body-centered cubic (BCC) oxygen sub-lattice but
differ in the dynamics of hydrogen atoms (protons)[2]. For the molecular
crystal ice-VII, the orientation of hydrogen-bonding is disordered and
continually changing as in hexagonal ice, obeying the ’ice rule’ [3]. When
temperature grows above thousands of Kelvin, ice-VII transforms into the
superionic phase VII” [4]. It has been suggested that the suitable conditions
for superionic ice lie deep inside the watery giants Uranus and Neptune and
may be common throughout the Universe [5, 1]. VII” is characterized by highly
mobile hydrogen ions (protons), behaving like a liquid and moving within the
BCC oxygen sub-lattice. The difference in the behavior of protons can result
in anomalies in thermodynamic and transport properties of ice.

The occurrence and geodynamic behaviors of these high-pressure ice polymorphs
(MPa-GPa range) have important effects on the thermal evolution of icy
planets. On this issue, thermal conductivity serves a key role for in-depth
understanding. However, despite the enormous phase transition regime and
proton transfer dynamics explored by previous efforts [6, 7, 8, 9, 1], it
seems likely that we still know little about the thermal transport properties
of dense ice across superionic transition.

Existing experimental efforts had been pursued to measure the thermal
conductivity of ice-VII up to 20 GPa [10, 11], but still far from the
condition of superionic regime due to the limitation of experimental
techniques under extreme conditions. From a theoretical point of view, the
dynamical nature of protons prevents the most commonly used tool, the lattice
dynamics approach, from tackling these issue. Another way to obtain thermal
conductivity is molecular dynamics simulation. However, ab initio method
requires an expensive computational cost to reach the long-time trajectories
with large simulation size required for estimation of correlation function
[12]. Moreover, the diverse local environments that characterize the different
relevant phases of water make classical force fields unfit for an accurate
simulation of their properties. Until now, the microscopic mechanism
determining the thermal conductivity of superionic ice at high pressure
remains unclear.

Figure 1:  (a) Atomic trajectories of ice-VII and VII” at different
temperatures during a 20-ps long run. The oxygen and hydrogen atoms are orange
and blue respectively. The cyan color is used to highlight the selected
hydrogen atoms that undergo transitions from bonded states between two
adjacent oxygen atoms to superionic states around different oxygen atoms. (b)
Temperature-dependent thermal conductivity $\kappa$ of ice-VII and VII” along
the isobar of P = 30 GPa. The gray vertical dashed line denotes the VII-VII”
phase boundary obtained from previous work [13].

Recent advances in machine-learning potential surface allow a full quantum-
mechanical, ab initio treatment of the interatomic interactions efficiently.
The deep potential water model is reported to predict a phase diagram close to
experiments [13], and its following applications have demonstrated its success
in estimating the thermal conductivity of water at extreme conditions [14,
15]. Therefore, in this work, we adopt the DP-SCAN water model and conduct a
series of deep potential molecular dynamics (DPMD) simulations to obtain the
thermal conductivity of VII and VII”, as well as diffusion coefficient,
spectral energy density (SED), and dynamic structure factor (DSF), to
understand the heat transport and to unravel the impact of mobile protons
across superionic phase transition.

Computational Details The DP model was trained with the DeePMD-kit package
[16, 17] using diverse ice crystal and liquid phase covering from ambient
condition to extreme thermodynamic state (p = 50 GPa, T = 2400 K). The
training data were obtained from density functional theory calculations using
the strongly constrained and appropriated normed (SCAN) exchange-correlation
functional. More details can be found in [13].

With DPMD simulation, the lattice thermal conductivity $\kappa$ is obtained
from the integration of the heat current autocorrelation function (HCACF),
known as Green–Kubo formula [18]. We performed a series of DPMD simulations
with the LAMMPS package [19]. A large supercell containing 1,296 atoms is used
to overcome the size effect (see Fig. S1 in the SI). The timestep was set to
0.5 fs and the Nosé-Hoover thermostat [20, 21] was employed in the NVT
ensemble. After a thermalization stage of 20 ps, the ensemble is switched into
the NVE ensemble to calculate the HCACF during the next 320 ps with the
correlation time set to 32 ps. To provide a representative sample for the
relevant statistical analysis, each (P, T) case repeats 20 times, with
independent initial velocity distribution. We note that such computation
complexity can hardly be achieved by the traditional AIMD method.

Non-monotonic behavior of thermal conductivity at elevated temperature. The
isobar of $p=30\ \rm{GPa}$ is chosen to investigate the thermal and proton
transport of ice-VII and VII”. As the temperature increases from 800 K to 1600
K, our DPMD simulations reproduce the superionic phase transition of ice
reported in previous works [7, 9, 13]. The atomic trajectories at temperatures
near the phase boundary are shown in Fig.1(a). More atomic trajectories can be
seen in Fig. S2 in SI. We can easily identify the different behaviors of
protons in the oxygen sub-lattice. At low temperatures, the hydrogen atoms are
bonded and only vibrate within the O-H···O bonds. At 1200 K close to the phase
boundary, the hydrogen atoms begin to initiate hopping between different
O-H···O bonds but remain bonded with oxygen atoms. The proton can migrate from
an O-H···O bond to another, leading to a fast change in the orientation of
water molecules. At higher temperatures, the protons diffuse freely out of the
bcc oxygen sub-lattice. Namely, the system transits into a superionic phase.

Correspondingly, $\kappa$ exhibits a non-monotonic trend, as shown in Fig.
1(b). Firstly, $\kappa$ decreases from $6.12\pm 0.14\ \rm{Wm^{-1}K^{-1}}$ at
800 K to a minimum value of $5.37\pm 0.21\ \rm{Wm^{-1}K^{-1}}$ at 1000 K. Then
a significant increase in $\kappa$ to a three-fold value of $17.28\pm 0.67\
\rm{Wm^{-1}K^{-1}}$ at T = 1200 K is observed. As the ice transits from VII
phase into VII” phase, $\kappa$ gradually increases and finally reaches a
plateau value of $22.72\pm 1.13\ \rm{Wm^{-1}K^{-1}}$ at 1600 K.

Figure 2: Temperature dependence of (a) thermal conductivity contributed by
heat conduction $\kappa_{conv}$ (b) diffusion coefficient along the isobar of
p = 30 GPa. Blue and orange markers denote the results of ice-VII and
superionic VII”, respectively. The gray vertical dashed line denotes the VII-
VII” phase boundary obtained from previous work [13].

Dominant role of proton diffusion in heat convection. The anomalous non-
monotonic trend of $\kappa$ is attributed to the significant change in proton
transport across the superionic transition. We extract the contributions of
heat convection and conduction to $\kappa$ ($\kappa_{conv}$ and
$\kappa_{cond}$) by decomposing the heat current into a heat convection term
and a heat conduction term, respectively (see more details in SI). As depicted
in Fig. 2(a), $\kappa_{conv}$ shows a monotonic increasing trend with
increasing temperature. At low temperatures, $\kappa_{conv}$ is close to zero.
At temperatures near the phase boundary, $\kappa_{conv}$ increases sharply. At
higher temperatures, $\kappa_{conv}$ gradually converges to a plateau similar
to $\kappa$. $\kappa_{conv}$ overwhelms the $\kappa$ after the superionic
transition, because of the fast diffusion of protons. The behavior of
$\kappa_{conv}$ dominates the dramatic increase and saturation of $\kappa$,
highlighting the importance of proton transfer dynamics in understanding the
heat transport in ice-VII and superionic VII”.

We plot the mean square displacements (MSD) of oxygen and hydrogen atoms in
Fig. S3 in SI. The diffusion coefficient can be estimated from the slope of
MSD from DPMD trajectories. Oxygen atoms have a flat curve of MSD and thus a
diffusion coefficient close to zero. In comparison, the MSD curves of proton
show diffusion characteristics and different slopes at different temperatures.
As shown in Fig. 2(b), the diffusion coefficient of proton $D_{H}$ increases
by two orders of magnitude, as the temperature increases from 1000 K to 1600
K. Besides, the diffusion behavior of proton at elevated temperatures can be
well described as an Arrhenius-like process, which is given as
$D(T)=Ae^{-E_{a}/k_{B}T}$, where $T$ is temperature, $A$ is a prefactor,
$E_{a}$ is the activation energy of the hopping mechanism leading to particle
diffusion, and $k_{B}$ is Boltzmann constant. By fitting the $D_{H}$ with the
Arrhenius model, the $E_{a}$ of ice-VII and VII” is $1.77\ \rm{eV}$ and $0.36\
\rm{eV}$ respectively. The much smaller $E_{a}$ of ice-VII” indicates the
weaker confinement and stronger diffusion of protons in superionic phase. The
temperature where the diffusion behavior of proton changes is consistent with
the phase boundary of superionic transition.

Figure 3:  Normalized dynamic structure factor $S(k,\omega)$ of hydrogen and
oxygen atoms along the isobar of p = 30 GPa at $k=0.07\ \rm{\mathbf{A}^{-1}}$.

Except for the mass diffusion, we further investigate the thermal diffusion
based on the dynamic structure factor $S(k,\omega)$ (see more details in SI).
The central Rayleigh peak encodes the thermal diffusion process as the
wavenumber is small enough to reach the hydrodynamic regime [22, 23]. At the
hydrodynamic limit, the shape of this central peak can be described by a
Lorentzian function with peak width relating to the thermal diffusivity
$D_{T}$, which gives $S(k,\omega)\propto
2D_{T}k^{2}/(\omega^{2}+(D_{T}k^{2})^{2})$.

We calculate the $S(k,\omega)$ of oxygen and hydrogen atoms. By fitting the
DPMD results with the hydrodynamic expression, we present the normalized
Rayleigh-Brillouin triplets for both H-contributed and O-contributed
$S(k,\omega)$ at a small wavenumber of $k=0.07\ \rm{\mathbf{A}^{-1}}$. As
presented in Fig. 3, as temperature increases, the broadening of central
Rayleigh peaks for protons is more significant than that for oxygen atoms,
indicating a more significant increase in thermal diffusion of protons
compared with that of oxygen atoms. The sharp increase in $\kappa$ across the
superionic transition stems from the dramatically enhanced diffusion of
protons and thus sharply increased thermal diffusion of protons.

Figure 4: Temperature dependence of (a) thermal conductivity contributed by
heat conduction $\kappa_{cond}$ and (b) normalized spectra energy density
$C(k,\omega)$, where the direction of wavevector $\mathbf{k}=(n_{x}\Delta
k_{x},0,0)$ is set to $x$ direction with the wavenumber resolution of $\Delta
k_{x}=2\pi/L_{x}\sim 0.07\ \rm{\mathring{A}^{-1}}$. The red and blue dotted
line denotes the linear dispersion relationship for longitudinal and
transverse acoustic branches, respectively.

Transverse mode softening across superionic transition. We also extract the
contributions of heat conduction to $\kappa$ ($\kappa_{cond}$). As presented
in Fig. 4(a), $\kappa_{cond}$ shows a monotonic decrease trend for ice-VII.
The $1/T$ dependence of $\kappa_{cond}$ is a typical characteristic that
reveals the dominant role of three-phonon scatterings [24]. At low
temperatures, $\kappa_{cond}$ is much larger than $\kappa_{conv}$, leading to
a decreasing trend. Therefore, the anomalous non-monotonic trend of $\kappa$
originates from the competing mechanism of heat conduction and convection.

We also note a discontinuity of $\kappa_{cond}$ near the phase boundary. On
one hand, it can be attributed to the density decrease accompanied by the VII-
VII” phase transition (see Fig. S5 in SI). On the other hand, a softening of
the transverse acoustic modes is observed. Here we calculate the spectra
energy density $C(k,\omega)$ from DPMD trajectories (see details in SI). The
$C(\mathbf{k},\omega)$ provides information on collective vibrational modes
(group velocity, lifetime) inside the complex ice polymorph. As shown in
Fig.4(b), we present the $C(k,\omega)$ in the low-frequency regime ($\nu\leq
20\ \rm{THz}$), which contains the longitudinal and transverse acoustic
branches that make dominant contributions to $\kappa_{cond}$. As the
wavenumber approaches the center of the first Brillouin zone, the dispersion
relationship exhibits a linear behavior, and the sound velocity can be
extracted. For ice-VII, the sound velocities of both longitudinal and
transverse acoustic branches decrease slightly with increasing temperature.
For ice-VII”, the sound velocity of longitudinal acoustic branches maintains a
slight decreasing trend while the sound velocity of transverse acoustic
branches shows a sudden larger decrease across the superionic transition. A
significant softening of the transverse acoustic modes across the superionic
transition is observed with decreased group velocity from $0.96\ \rm{km/s}$ to
$0.74\ \rm{km/s}$. Moreover, as temperature increases, the vibrational energy
peaks of both longitudinal and transverse acoustic branches exhibit
broadening, corresponding to the reduction in phonon lifetimes.

Overall, by summarizing our results, we can understand the anomalous behavior
of $\kappa$ in ice-VII across the superionic transition. At moderate
temperatures, the propagation of lattice vibration mode dominates the $\kappa$
and exhibits a typical $1/T$ dependence due to the three-phonon scattering
process. At T = 1000 K, the onset of hydrogen diffusion between O-H···O pairs
leads to an exponential increase by two orders of magnitude in the diffusion
coefficient of protons. These protons hop within the oxygen-formed BCC sub-
lattice and carry heat, creating a non-negligible contribution via heat
convection. As temperature increases above the superionic transition threshold
(T = 1250 K), the first-order phase transition occurs. The significantly
enhanced mobility of hydrogen, combined with softening in transverse acoustic
branches, leads to a saturated value above T = 1300 K. Under such conditions,
the diffusing protons can experience a stronger thermal diffusion process as
compared with oxygen atoms.

In conclusion, by utilizing the newly developed DP-SCAN water model, we
investigate the microscopic mechanism behind the anomalous thermal and proton
transport of high-pressure ice across the superionic transition. We explain
the anomalous increasing trend of $\kappa$ with elevated temperature and
illustrate the important role of proton diffusion in superionic ice. To
overcome the limitation of the traditional lattice dynamics approach which
requires high-order force constants to correct the quasi-harmonic
approximation, here we extract all the mass transport, heat transport, and
collective dynamics from long-time large-scale molecular dynamics
trajectories, without any assumptions or approximations made. These approaches
combined with ab initio accurate deep neural network potential energy surface
model, can be applied to various complex materials, including proton-disorder
ice polymorphs, proton-diffusion superionic crystal, and amorphous materials.

## Acknowledgment

This work was supported by the National Key R&D Program of China under Grant
No. 2017YFA0403200, the NSAF under Grant No. U1830206, the Science and
Technology Innovation Program of Hunan Province under Grant No. 2021RC4026.

## References

  * Prakapenka  _et al._ [2021] V. B. Prakapenka, N. Holtgrewe, S. S. Lobanov, and A. F. Goncharov, Structure and properties of two superionic ice phases, Nature Physics 17, 1233–1238 (2021).
  * Kang  _et al._ [2013] D. Kang, J. Dai, H. Sun, Y. Hou, and J. Yuan, Quantum simulation of thermally-driven phase transition and oxygen k-edge x-ray absorption of high-pressure ice, Sci Rep 3, 3272 (2013).
  * Bernal and Fowler [1933] J. D. Bernal and R. H. Fowler, A theory of water and ionic solution, with particular reference to hydrogen and hydroxyl ions, J Chem Phys 1, 515 (1933).
  * Cavazzoni  _et al._ [1999] C. Cavazzoni, G. L. Chiarotti, S. Scandolo, E. Tosatti, M. Bernasconi, and M. Parrinello, Superionic and metallic states of water and ammonia at giant planet conditions, Science 283, 44 (1999).
  * Chaplin [2008] M. Chaplin, Water structure and science, (2008).
  * Schwegler  _et al._ [2008] E. Schwegler, M. Sharma, F. Gygi, and G. Galli, Melting of ice under pressure, Proc Natl Acad Sci U S A 105, 14779 (2008).
  * Hernandez and Caracas [2016] J.-A. Hernandez and R. Caracas, Superionic-superionic phase transitions in body-centered cubic ${\mathrm{h}}_{2}\mathrm{O}$ ice, Phys. Rev. Lett. 117, 135503 (2016).
  * Hernandez and Caracas [2018] J. A. Hernandez and R. Caracas, Proton dynamics and the phase diagram of dense water ice, J Chem Phys 148, 214501 (2018).
  * Queyroux  _et al._ [2020] J.-A. Queyroux, J.-A. Hernandez, G. Weck, S. Ninet, T. Plisson, S. Klotz, G. Garbarino, N. Guignot, M. Mezouar, M. Hanfland, J.-P. Itié, and F. Datchi, Melting curve and isostructural solid transition in superionic ice, Phys. Rev. Lett. 125, 195501 (2020).
  * Andersson and Inaba [2005] O. Andersson and A. Inaba, Thermal conductivity of crystalline and amorphous ices and its implications on amorphization and glassy water, Phys Chem Chem Phys 7, 1441 (2005).
  * Chen  _et al._ [2011] B. Chen, W.-P. Hsieh, D. G. Cahill, D. R. Trinkle, and J. Li, Thermal conductivity of compressed h2o to 22 gpa: A test of the leibfried-schlömann equation, Physical Review B 83, 10.1103/PhysRevB.83.132301 (2011).
  * Grasselli  _et al._ [2020] F. Grasselli, L. Stixrude, and S. Baroni, Heat and charge transport in h2o at ice-giant conditions from ab initio molecular dynamics simulations, Nat Commun 11, 3605 (2020).
  * Zhang  _et al._ [2021] L. Zhang, H. Wang, R. Car, and W. E, Phase diagram of a deep potential water model, Phys Rev Lett 126, 236001 (2021).
  * Zhang  _et al._ [2023] C. Zhang, M. Puligheddu, L. Zhang, R. Car, and G. Galli, Thermal conductivity of water at extreme conditions, J Phys Chem B 127, 7011 (2023).
  * Yang  _et al._ [2022] F. Yang, Q. Zeng, B. Chen, D. Kang, S. Zhang, J. Wu, X. Yu, and J. Dai, Lattice thermal conductivity of mgsio3 perovskite and post-perovskite under lower mantle conditions calculated by deep potential molecular dynamics, Chinese Physics Letters 39, 116301 (2022).
  * Wang  _et al._ [2018] H. Wang, L. Zhang, J. Han, and E. Weinan, Deepmd-kit: A deep learning package for many-body potential energy representation and molecular dynamics, Computer Physics Communications 228, 178 (2018).
  * Zeng  _et al._ [2023] J. Zeng, D. Zhang, D. Lu, P. Mo, Z. Li, Y. Chen, M. Rynik, L. Huang, Z. Li, S. Shi, Y. Wang, H. Ye, P. Tuo, J. Yang, Y. Ding, Y. Li, D. Tisi, Q. Zeng, H. Bao, Y. Xia, J. Huang, K. Muraoka, Y. Wang, J. Chang, F. Yuan, S. L. Bore, C. Cai, Y. Lin, B. Wang, J. Xu, J. X. Zhu, C. Luo, Y. Zhang, R. E. A. Goodall, W. Liang, A. K. Singh, S. Yao, J. Zhang, R. Wentzcovitch, J. Han, J. Liu, W. Jia, D. M. York, W. E, R. Car, L. Zhang, and H. Wang, Deepmd-kit v2: A software package for deep potential models, J Chem Phys 159, 054801 (2023).
  * Mcquarrie [1965] D. Mcquarrie, Statistical mechanics (1965).
  * Plimpton [1995] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, Journal of computational physics 117, 1 (1995).
  * Nosé [1984] S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, The Journal of chemical physics 81, 511 (1984).
  * Hoover [1985] W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A 31, 1695 (1985).
  * Hansen and McDonald [1990] J.-P. Hansen and I. R. McDonald,  _Theory of simple liquids_ (Elsevier, 1990).
  * Zeng  _et al._ [2021] Q. Zeng, X. Yu, Y. Yao, T. Gao, B. Chen, S. Zhang, D. Kang, H. Wang, and J. Dai, Ab initio validation on the connection between atomistic and hydrodynamic description to unravel the ion dynamics of warm dense matter, Phys. Rev. Research 3, 033116 (2021).
  * Roufosse and Klemens [1973] M. Roufosse and P. G. Klemens, Thermal conductivity of complex dielectric crystals, Physical Review B 7, 5379 (1973).


11institutetext: Indian Institute of Technology Guwahati, Guwahati, India
11email: {saswatajana<EMAIL_ADDRESS>22institutetext: Luiss University,
Rome, Italy
22email<EMAIL_ADDRESS>22email<EMAIL_ADDRESS>22email:
<EMAIL_ADDRESS>University of Ioannina, Ioannina,
Greece
33email<EMAIL_ADDRESS>
# Online Drone Scheduling for Last-mile Delivery

Saswata Jana 11 Giuseppe F. Italiano 22 Manas Jyoti Kashyop 22 Athanasios L.
Konstantinidis 22 Evangelos Kosinas 33 Partha Sarathi Mandal 11

###### Abstract

Delivering a parcel from the nearest distribution hub to the customer’s
doorstep is called the last-mile delivery step in delivery logistics. Last-
mile delivery is the costliest and most time-consuming phase of the entire
delivery logistics. It is gaining importance due to the substantial rise in
the e-commerce market and customers’ expectations of quicker delivery.
Delivery by unmanned aerial vehicles (commonly known as drones) has garnered
significant interest due to its environment-friendly operations, faster
delivery times, low operational charges, and the ability to reach remote
areas. In this paper, we study a hybrid truck-drones model for the last-mile
delivery step. A truck moves on a path carrying drones and parcels. Drones do
the deliveries, whereas the moving truck is used as a base for launching and
landing drones synchronously. We define the online drone scheduling problem,
where the truck moves in a predefined path, and the customer’s requests appear
online during the truck’s movement. The objective is to schedule a drone
associated with every request to minimize the number of drones used subject to
the battery budget of the drones and compatibility of the schedules. We
propose a 3-competitive deterministic algorithm using the next-fit strategy
and 2.7-competitive algorithms using the first-fit strategy for the problem
with $O(\log n)$ worst-case time complexity per request, where $n$ is the
maximum number of active requests at any time. We also introduce online
variable-size drone scheduling problem (OVDS). Here, we know all the
customer’s requests in advance; however, the drones with different battery
capacities appear online. The objective is to schedule customers’ requests for
drones to minimize the number of drones used. We propose a
$(2\alpha+1)$-competitive algorithm for the OVDS problem with total running
time $O(n\log n)$ for $n$ customer requests, where $\alpha$ is the ratio of
the maximum battery capacity to the minimum battery capacity of the drones.
Finally, we address how to generate intervals corresponding to each customer
request when there are discrete stopping points on the truck’s route, from
where the drone can fly and meet with the truck.

###### Keywords:

Online Algorithm Optimization Drone-Delivery Scheduling Last-mile Delivery

##  1 Introduction

Rapid developments in unmanned aerial vehicle (UAV) technology, popularly
known as drone technology, motivate scientists and researchers to use drones
in various operations. The flexibility and mobility of drones allow them to be
used in healthcare, weather forecasts, surveillance, disaster management,
agriculture, artificial intelligence, and many more. Recently, delivery giant
companies, e.g., Amazon, DHL, and FedEx, have started using drones for the
last-mile delivery [2], [8]. Last-mile delivery is the delivery logistic
system’s final and costliest step. In this step, the package is collected from
the distribution hub and delivered to the customer using a vehicle, e.g., a
delivery truck. The problem of delivering goods using vehicles is studied as
the routing problem, where the objective is to minimize the total make-span.
If a single vehicle without any capacity constraint is used in the routing
problem, the problem reduces to the traveling salesman problem (TSP). If
several vehicles with capacity constraints are used, the routing problem is
called the vehicle routing problem (VRP). In literature, several variants of
both problems are considered, e.g., deployment of depots for refilling the
delivery trucks, bounding the number of depots, and introducing time windows
to serve a customer [19], [20], all with the final objective to minimize the
total make-span. The use of drones along with delivery vehicles was studied
for the first time by Murray and Chu [21], where the authors proposed a model
called the flying sidekicks traveling salesman problem (FSTSP). In this model,
a delivery truck and a single drone cooperate to deliver the packages. The
truck starts its journey from the depot along with the deliverables and
carries the drone. A customer is either served by the truck or by the drone.
The authors presented a mixed integer linear programming (MILP) formulation
and heuristics for finding the best schedule using the model. Later, the
authors in [6], [22] modified the FSTSP model and introduced multiple drones
(extended FSTSP model). However, in the extended FSTSP model, the authors did
not consider the limitation of the flight endurance of the drones. Daknama and
Kraus [7] proposed a variation of the extended FSTSP model where every drone
has a limited battery capacity (budgeted extended FSTSP model) and thus
returns to the truck after completing the delivery for recharging. Authors in
[24] considered the budgeted extended FSTSP model from the worst-case point of
view and gave comparative results depending on the number of drones and the
relative speed of drone and truck. In [1], [12],[17], and [23], authors
considered different UAV-based frameworks with a focus on battery recharging
policy. They proposed heuristics-based approaches to optimize different
objective functions, such as minimizing the total travel cost, minimizing the
energy consumption, and maximizing the total revenue or coverage efficiency.

The use of drones in delivery logistics has several advantages. Reducing
$CO_{2}$ emissions is a major concern in the transportation system. Recent
studies suggest that the use of drones in the delivery logistic system reduces
$CO_{2}$ emissions compared to the use of only delivery trucks or human-
operated ground vehicles, for example: [5], [13], and [18]. Using drones in
the delivery logistics system has other motivations, such as avoiding traffic
congestion by flying, following a shorter route toward the customer’s
doorstep, and contactless delivery. Delivery with minimal human involvement or
contactless delivery became crucial after the recent pandemic. However, using
only drones has certain limitations, such as the size of the deliverable
packages and the distance of the customer from the launching location due to
limited battery capacity. Therefore, a hybrid model with a delivery truck and
several drones (truck-drones model) is more desirable [24]. In the truck-
drones model, the truck containing all the drones and the packages follows
conventional ground routes, commencing its journey from the warehouse. A
customer will be delivered either by truck or by drone. In [4], authors
considered that the truck has some specific discrete stopping points where the
drone can take off from the truck and land on the truck. A customer places a
request with information about its location, and based on that, a launching
point and a rendezvous point for the drone to be used to serve the request are
computed. Several works (e.g., [14]) considered that customer requests will
appear online, and the objective is to minimize the total make-span. Sorbeli
et al. [3] proposed a similar model in the offline version, where all the
customer positions and the corresponding launching and meeting points are
known in advance. The objective is to maximize the profit by serving a subset
of customers with a fixed number of drones subject to battery and
compatibility constraints. Authors in [15] took a variant of the above model
where the truck has a sufficient number of drones to complete all the
deliveries, and the objective is to use the minimum number of drones. The
authors in [3], [15], and [16] proposed several heuristics and approximation
algorithms for this offline version.

In this work, we follow the hybrid approach and consider the truck-drones
model. Our objective is to use the minimum number of drones. For simplicity,
we do not mention the customer requests that are close to the truck path and
are served by the truck. Further, we consider that the truck has some specific
discrete stopping points where the drone can take off from the truck and land
on the truck. We compute the launch and rendezvous point of each customer
request for a drone, and these requests appear online. Based on the launch and
rendezvous point, an interval generator generates an interval for every
customer request. We call the amount of battery consumption of a drone used to
serve a specific request as the cost of the interval corresponding to that
request. The objective is to use the minimum number of drones to serve all the
requests. Thus, we differ from the previous works that considered online
customer requests intending to minimize the total make-span. We assumed that
the truck moves on a known path containing a sufficient number of drones, and
the path has some discrete stops from where the drone can launch or meet the
truck. Also, the time to reach or leave those points by truck is known. Now a
request from a customer will come at any point on the path. We need to select
the best points for launching and rendezvous so that the cost to complete the
delivery is minimum.

###  1.1 Our Contribution

  * •

We define the online drone scheduling problem to minimize the number of drones
in the last-mile delivery scenario.

  * •

We propose a deterministic 3-competitive algorithm for the online drone
scheduling problem by using the next-fit bin packing strategy.

  * •

We propose a 2.7-competitive algorithm for the problem by applying the first-
fit bin packing strategy.

  * •

We show that all the above algorithms use $O(\log n)$ worst-case time
complexity per scheduling, where $n$ is the maximum number of active requests
at any time. We define some data structures to achieve this time efficiency.

  * •

We define another variant of the problem named online variable-size drone
scheduling and propose an $(2\alpha+1)$\- competitive algorithm with time
complexity $O(n\log n)$, where $\alpha$ is the ratio between the maximum and
minimum battery capacity of the drones.

  * •

We also discuss how to generate delivery intervals in an online setting having
some discrete stop points on the truck path.

###  1.2 Preliminaries

A graph is an _interval graph_ if there is a bijection between its vertices
and a family of closed intervals of the real line such that two vertices are
adjacent if and only if the two corresponding intervals intersect. Such a
bijection is called an _interval representation_ of the graph, denoted by
$\mathcal{I}$. We identify the intervals of the given representation with the
vertices of the graph, interchanging these notions appropriately. Whether a
given graph is an interval graph can be decided in linear time and if so, an
interval representation can be generated in linear time [10]. Notice that
every induced subgraph of an interval graph is an interval graph. A clique of
$G$ is a set of pairwise adjacent vertices of $G$, and a maximal clique of $G$
is a clique of $G$ that is not properly contained in any clique of $G$. An
independent set of $G$ is a set of pairwise non-adjacent vertices of $G$.

Organization: The rest of the paper is organized as follows. Section 2
introduces the problem definition for the online drone scheduling problem and
a few algorithms to solve it. Section 3 defines the online variable-size drone
scheduling problem along with an algorithm. Section 4 discusses interval
generation. Finally we conclude in Section 5.

##  2 Drone scheduling through bin packing with conflicts

We consider the interval representation of the underlying interval graph $G$.
Let $I_{t}=[l_{t},r_{t}]$ be the interval inserted at time $t$. Let $v_{t}$ be
the vertex in $G$ corresponding to $I_{t}$. Therefore, the insertion or
deletion of an interval is the insertion or deletion of the corresponding
vertex of the underlying graph $G$. In our work, every interval $I_{t}$ is
also associated with a cost (${\sf cost}(I_{t})$) that corresponds to the
amount of battery consumption of a drone while serving the customer request
corresponding to the interval. In this section, we consider the online drone
scheduling problem as proper coloring of intervals with a fixed budget $B$ for
every color, where $B$ is a positive integer and is part of the input. A color
assigned to an interval corresponds to a drone assigned to serve the request
represented by the interval and budget $B$ of the color corresponds to the
battery capacity of the drone.

We define the following problem.

###### Definition 1

(Online Drone Scheduling) Intervals are appearing online. Suppose at time step
$t$, interval $I_{t}=[l_{t},r_{t}]$ appears with cost ${\sf cost}(I_{t})$ with
$l_{t}\geq t$. The goal is to assign a color to the interval with the
following constraints:
$1$. Any two Intervals with the same color are disjoint (non-overlapping).
$2$. For a particular color, the summation of the costs of all intervals with
that color is at most $B$.

To solve the problem we will use the algorithms developed for online bin
packing with conflicts in [9]. For interval $I_{t}$ at time $t$, our goal is
to compute a color ${\sf Color}(I_{t})$. A color ${\sf Color}(I_{t})$ is a
$2$-tuple, ${\sf Color}(I_{t})$ = $({\sf idNumber}(I_{t}),{\sf
binNumber}(I_{t}))$, with the following properties.

  1. 1.

(conflict free) No two overlapping intervals have the same ${\sf idNumber}$.
Thus, the ${\sf idNumber}$ ensures a conflict free property.

  2. 2.

(budget constraint) For all the intervals with the same ${\sf idNumber}$ and
${\sf binNumber}$ (in other words the same color), the summation of their
costs must be at most $B$. Therefore, for the intervals with the same ${\sf
idNumber}$, ${\sf binNumber}$ ensures that the budget constraint associated
with every color is satisfied.

  3. 3.

(proper coloring) The above two properties imply that for any two intervals
$I_{t}\neq I_{t^{\prime}}$ and $I_{t}\cap I_{t^{\prime}}\neq\emptyset$, $({\sf
idNumber}(I_{t}),{\sf binNumber}(I_{t}))\neq({\sf idNumber}(I_{t^{\prime}}),$
${\sf binNumber}(I_{t^{\prime}}))$.

Based on the above properties, we use the following to compute the color.

  1. 1.

${\sf idNumber}$ is computed using an online interval coloring algorithm where
the colors are represented using positive integer numbers. Note that the
coloring algorithm used at this step does not have any budget associated with
the colors and thus can be assigned to any number of non-overlapping
intervals. For an interval $I_{t}$, the color (the non-negative integer)
allotted to $I_{t}$ at this step using the online coloring algorithms is
assigned as the ${\sf idNumber}(I_{t})$.

  2. 2.

Once the ${\sf idNumber}(I_{t})$ is computed for an interval $I_{t}$ in the
previous step, next step is to compute ${\sf binNumber}(I_{t})$ ensuring that
all the intervals with the same color as $({\sf idNumber}(I_{t}),{\sf
binNumber}(I_{t}))$ satisfies the budget constraint $B$. Therefore, we use an
online bin packing algorithm to compute the ${\sf binNumber}$ for $I_{t}$
considering all the previously inserted intervals with the ${\sf idNumber}$
equal to ${\sf idNumber}(I_{t})$.

For an interval $I_{t^{\prime}}=[l_{t^{\prime}},r_{t^{\prime}}]$, color of the
interval is computed at time $t=l_{t^{\prime}}$. Further, we assume that
intervals have distinct endpoints. Therefore, at the time step $t$, the
following scenarios occur while handling the update.

  1. 1.

If a new interval $I_{t}$ appears, then $I_{t}$ is inserted. Further, at time
step $t$ at most one interval appears.

  2. 2.

If there is an existing interval $I_{t^{\prime}}$ with $t^{\prime}<t$ and
$r_{t^{\prime}}=t$, then $I_{t^{\prime}}$ is deleted.

  3. 3.

If there is an interval $I_{t^{\prime}}$ with $t^{\prime}\leq t$ and
$l_{t^{\prime}}=t$, then color for $I_{t^{\prime}}$ is computed.

  4. 4.

As intervals have distinct endpoints, scenario $2$ and $3$ are mutually
exclusive.

Because of the properties of our original online drone scheduling problem, the
intervals appearing online satisfy the following properties.

  * •

Interval appearing at time step $t$ has left endpoint at least $t$.

  * •

Our deletion scheme is that at time step $t$, if there is an interval
$I_{t^{\prime\prime}}$ with $t^{\prime\prime}<t$ and $r_{t^{\prime\prime}}=t$,
then we delete $I_{t^{\prime\prime}}$. Otherwise, if there is an interval
$I_{t^{\prime}}$ with $t^{\prime}\leq t$ and $l_{t^{\prime}}=t$, the color for
$I_{t^{\prime}}$ is computed. Note that because of our deletion scheme, while
computing the color for $I_{t^{\prime}}$, all the intervals overlapping with
$I_{t^{\prime}}$ must contain the point $t$ forming a clique. Further, because
of our coloring scheme, colors are assigned to the intervals in the increasing
order of their left endpoints.

###  2.1 Data structures

We use the following data structures.

  * •

We store every interval $I_{t}=[l_{t},r_{t}]$ with ${\sf cost}(I_{t})$ in a
map $\mathcal{I}$ using $t$ as key.

  * •

The left endpoint $l_{t}$ for every interval $I_{t}=[l_{t},r_{t}]$ is stored
as $(l_{t},t)$ in a min heap $\sf H^{min}_{L}$ using $l_{t}$ as key.

  * •

The right endpoint $r_{t}$ for every interval $I_{t}=[l_{t},r_{t}]$ is stored
as $(r_{t},t)$ in a min heap $\sf H^{min}_{R}$ using $r_{t}$ as key.

  * •

Color of every interval $I_{t}=[l_{t},r_{t}]$ is stored as $({\sf
idNumber}(I_{t}),{\sf binNumber}(I_{t}))$ in a map ${\sf Color}$ with $t$ as
the key. In the case, interval $I_{t}$ is deleted or the color for $I_{t}$ is
yet to be computed, ${\sf Color}[t]$ stores NULL.

The operation $\sf extractMin()$ returns the minimum element in a min heap.
The operation $\sf insertHeap()$ inserts an element into a heap and $\sf
deleteHeap()$ deletes the root element of the heap (that is the minimum
element in the case of min heap and the maximum element in the case of max
heap). $\sf searchMap()$ is used to search for a key in a map and $\sf
insertMap()$ is used to insert an entry into a map. In addition, we will use
some data structures specific to the computation of ${\sf idNumber}$ and ${\sf
binNumber}$, and we will describe them in the relevant sections.

###  2.2 Computing ${\sf idNumber}$

At time step $t$, we consider the interval
$I_{t^{\prime}}=[l_{t^{\prime}},r_{t^{\prime}}]$ with $l_{t^{\prime}}=t$ and
$t^{\prime}\leq t$ to compute ${\sf idNumber}(I_{t^{\prime}})$. Let
$\mathcal{I}_{t}$ be the set of intervals overlapping with $I_{t^{\prime}}$
and contain the point $t$. Since all the intervals have distinct endpoints,
the left endpoint of every interval in $\mathcal{I}_{t}$ is strictly smaller
than $t$. We observe the following properties about the set of intervals
$\mathcal{I}_{t}$.

  1. 1.

$\mathcal{I}_{t}$ forms a clique because all the intervals in
$\mathcal{I}_{t}$ contains the point $t$.

  2. 2.

Because of our deletion scheme, $\mathcal{I}_{t}$ is the only set of intervals
overlapping with the interval $I_{t^{\prime}}$ and with left endpoints
strictly less than $t$.

  3. 3.

As of our coloring scheme, color for each interval in $\mathcal{I}_{t}$ is
already computed.

  4. 4.

Any interval not in $\mathcal{I}_{t}$ and overlapping with $I_{t^{\prime}}$
must have left endpoint strictly greater than $t$. Because of our coloring
scheme, the color for such an interval is not yet computed and we do not
consider such an interval while computing the color for $I_{t^{\prime}}$ in
the current step.

Our goal is to find the smallest ${\sf idNumber}$ not assigned to any interval
in $\mathcal{I}_{t}$ and assign it as ${\sf idNumber}(I_{t^{\prime}})$. Along
with the data structures described in Section 2.1, we use the following data
structures specifically to compute ${\sf idNumber}$.

  * •

A counter ${\sf gIdNumber}$ that stores the maximum value of ${\sf idNumber}$
used so far. It is initialized to $0$.

  * •

A min heap $\sf H^{min}_{F{\sf idNumber}}$ (a heap of free ${\sf idNumber}$s)
containing the numbers in the range $\\{1,\cdots,{\sf gIdNumber}\\}$ that are
not used by any existing interval as ${\sf idNumber}$.

  * •

A balanced binary search tree $\sf T_{U{\sf idNumber}}$ containing the numbers
in the range $\\{1,\cdots,{\sf gIdNumber}\\}$ that are used by an existing
interval as ${\sf idNumber}$. The operation $\sf insertBST()$ inserts an
element to the tree and $\sf deleteBST()$ deletes an element from the tree.

Note that because of our deletion scheme, the only ${\sf idNumber}$s in use
are the ones corresponding to the colors of the intervals in
$\mathcal{I}_{t}$. Therefore, the tree $\sf T_{U{\sf idNumber}}$ stores only
the ${\sf idNumber}$s of the intervals in $\mathcal{I}_{t}$ maintaining the
exact set of ${\sf idNumber}$s that are forbidden for us to use. Similarly,
the heap $\sf H^{min}_{F{\sf idNumber}}$ maintains all the free ${\sf
idNumber}$s that are potential candidates for us to use. Using the above data
structures, we design the procedure Algorithm 1.

1 if  _$\sf H^{min}_{F{\sf idNumber}}$ is Empty_ then

2 ${\sf gIdNumber}\leftarrow{\sf gIdNumber}+1$

3 ${\sf idNumber}(I_{t^{\prime}})\leftarrow{\sf gIdNumber}$

4

5else

6 ${\sf idNumber}(I_{t^{\prime}})\leftarrow\sf extractMin(\sf H^{min}_{F{\sf
idNumber}})$

7 $\sf deleteHeap(\sf H^{min}_{F{\sf idNumber}})$

8

9$\sf insertBST(\sf T_{U{\sf idNumber}},{\sf idNumber}(I_{t^{\prime}}))$

Algorithm 1 compute-idNumber$(I_{t^{\prime}})$

###### Lemma 1

compute-idNumber() compute ${\sf idNumber}$ for each interval
$I_{t^{\prime}}$.

###### Proof

Procedure compute-idNumber() (refer to Algorithm 1), first checks if there is
a previously used ${\sf idNumber}$ available for reuse. In the case of such
availability, it selects the minimum one and assigns it as the ${\sf
idNumber}$ of the current interval, and adjusts the relevant data structures
suitably (Line $4-6$ in Algorithm 1). Otherwise, all the intervals overlapping
with the current interval and excluding the current interval form a clique of
size ${\sf gIdNumber}$. Therefore, a new ${\sf idNumber}$ must be introduced
to assign it as the ${\sf idNumber}$ of the current interval (Lines $1-3$ in
Algorithm 1). Therefore, the procedure compute-idNumber() always assigns to
the current interval the least ${\sf idNumber}$ not used by any of the
existing intervals overlapping with it.

###### Lemma 2

Running time of procedure compute-idNumber() is $O(\log n)$ in the worst-case.

###### Proof

The running time of procedure compute-idNumber() is dominated by the standard
operations in min heap and binary search tree. Therefore, in the worst-case
the procedure takes $O(\log n)$ time, where $n$ is the maximum number of
elements present in the heap or the binary search tree throughout the handling
of the online interval sequence.

###### Lemma 3

${\sf gIdNumber}\leq\sf OPT_{off}$, where $\sf OPT_{off}$ is the number of
colors used by an optimal offline interval coloring algorithm without any
budget constraints on the colors.

###### Proof

The optimal offline interval coloring uses a greedy strategy. It considers the
intervals in the sorted order of their left endpoints and assigns an interval
to the least numbered color not already assigned to any of its neighbors.
Notice that the procedure compute-idNumber() mimics the offline greedy
strategy. Further, because of our deletion scheme, the set of free ${\sf
idNumber}$s available is a super set of the set of available colors in the
offline greedy strategy. Since ${\sf gIdNumber}$ denotes the maximum value of
${\sf idNumber}$ used so far by the algorithm, we have ${\sf gIdNumber}\leq\sf
OPT_{off}$.

###  2.3 Computing ${\sf binNumber}$

At time step $t$, after computing ${\sf idNumber}(I_{t^{\prime}})$ for the
interval $I_{t^{\prime}}=[l_{t^{\prime}},r_{t^{\prime}}]$ with
$l_{t^{\prime}}=t$ and $t^{\prime}\leq t$, the next task is to compute ${\sf
binNumber}$ for $I_{t^{\prime}}$. First we use the next-fit strategy from
online bin packing to compute the ${\sf binNumber}$.

###### Definition 2

(next-fit) The next-fit algorithm for online bin packing maintains an active
bin. Initially, the active bin is empty. When a new item appears, if the item
fits into the active bin, then the next-fit algorithm places the item in the
active bin. Otherwise, the current active bin is closed, a new bin is marked
as active, and the current item is placed in the new active bin.

To use the next-fit strategy for our purpose, we maintain an active bin number
corresponding to every ${\sf idNumber}$ along with its capacity (precisely,
remaining capacity). Once an item with a specific ${\sf idNumber}$ does not
fit into the current active bin number, a new bin is made active with capacity
$B$. We use the following data structures to compute ${\sf binNumber}$ using
next-fit.

  * •

A map $\sf idTobinNumber$ with ${\sf idNumber}$ as key which stores the
currently active bin number corresponding to the ${\sf idNumber}$. For every
new ${\sf idNumber}$, the entry $\sf idTobinNumber[{\sf idNumber}]$ is
initialized to $1$ (that is the bins corresponding to a particular ${\sf
idNumber}$ are numbered starting from $1$).

  * •

A map $\sf binToCapacity$ with ${\sf idNumber}$ as the key which stores the
remaining capacity of the currently active bin number corresponding to the
${\sf idNumber}$. For every new ${\sf idNumber}$, the entry $\sf
binToCapacity[{\sf idNumber}]$ is set to $B$.

At time step $t$, using the above data structures, we design the following
procedure (Algorithm 2) to compute ${\sf binNumber}$ for an interval
$I_{t^{\prime}}=[l_{t^{\prime}},r_{t^{\prime}}]$ with $t^{\prime}\leq t$ and
$l_{t^{\prime}}=t$ after computing ${\sf idNumber}(I_{t^{\prime}})$.

1 if  _$\sf searchMap(\sf idTobinNumber,{\sf idNumber}(I_{t^{\prime}}))$ is
Empty_ then

2 $\sf insertMap(\sf idTobinNumber,<{\sf idNumber}(I_{t^{\prime}}),1>)$

3 $\sf insertMap(\sf binToCapacity,<{\sf idNumber}(I_{t^{\prime}}),B-{\sf
cost}(I_{t^{\prime}})>)$

4

5else

6 if  _$\sf binToCapacity[{\sf idNumber}(I_{t^{\prime}})]\geq{\sf
cost}(I_{t^{\prime}})$_ then

7 $\sf binToCapacity[{\sf idNumber}(I_{t^{\prime}})]\leftarrow\sf
binToCapacity[{\sf idNumber}(I_{t^{\prime}})]-{\sf cost}(I_{t^{\prime}})$

8

9 else

10 $\sf idTobinNumber[{\sf idNumber}(I_{t^{\prime}})]\leftarrow\sf
idTobinNumber[{\sf idNumber}(I_{t^{\prime}})]+1$

11 $\sf binToCapacity[{\sf idNumber}(I_{t^{\prime}})]\leftarrow B-{\sf
cost}(I_{t^{\prime}})$

12

13${\sf binNumber}(I_{t^{\prime}})\leftarrow\sf idTobinNumber[{\sf
idNumber}(I_{t^{\prime}})]$

Algorithm 2 compute-binNumber-Nextfit$(I_{t^{\prime}})$

###### Lemma 4

compute-binNumber-Nextfit() computes ${\sf binNumber}$ for each
$I_{t^{\prime}}$.

###### Proof

For the interval $I_{t^{\prime}}$, the procedure compute-binNumber-Nextfit()
(refer to Algorithm 2) first checks if there is any active bin corresponding
to ${\sf idNumber}(I_{t^{\prime}})$. In the case when there is no active bin
for ${\sf idNumber}(I_{t^{\prime}})$, the first active bin with number $1$ is
created, the remaining capacity of the new bin is adjusted suitably, and ${\sf
binNumber}(I_{t^{\prime}})$ is assigned the value $1$ (Line $1-3$ and Line
$10$). Otherwise, if the current active bin has sufficient remaining capacity
to accommodate interval $I_{t^{\prime}}$, then the remaining capacity of the
bin is adjusted suitably, and ${\sf binNumber}(I_{t^{\prime}})$ is assigned
the number of the current active bin (Line $5-6$ and Line $10$). In the case
when the remaining capacity of the current active bin is insufficient to
accommodate $I_{t^{\prime}}$, a new bin corresponding to ${\sf
idNumber}(I_{t^{\prime}})$ is created, and is assigned as ${\sf
binNumber}(I_{t^{\prime}})$ after appropriately adjusting the remaining
capacity (Line $7-9$ and Line $10$). Therefore, procedure compute-binNumber-
Nextfit() mimics the next-fit strategy ensuring the constrain that sum of the
costs of all intervals with color $({\sf idNumber}(I_{t^{\prime}}),{\sf
binNumber}(I_{t^{\prime}}))$ is at most $B$.

###### Lemma 5

The running time of procedure compute-binNumber-Nextfit() is $O(\log n)$ in
the worst-case.

###### Proof

The running time of procedure compute-binNumber-Nextfit() (Algorithm 2) is
dominated by the standard search ($\sf searchMap()$) and insert ($\sf
insertMap()$) operation in a map data structure. The operation $\sf
searchMap()$ takes $O(1)$ time in the worst-case. The operation $\sf
insertMap()$ takes $O(\log n)$ time in the worst-case where $n$ is the maximum
number of elements present in the map data structure throughout handling the
online interval sequence. Therefore, the worst-case running time of procedure
compute-binNumber-Nextfit() is $O(\log n)$.

Finally, we present the algorithm $\sf DSusingBPC$ (online drone scheduling
using online bin packing with conflicts, refer to Algorithm 3 for pseudo code)
which handles the update at time $t$ as follows.

  * •

If a new interval $I_{t}$ is to be inserted, then insertion is performed by
updating the relevant map and heap data structures (Line $1-4$ in Algorithm
3).

  * •

If there is an interval $I_{t^{\prime\prime}}$ with $r_{t^{\prime\prime}}=t$,
then $I_{t^{\prime\prime}}$ is deleted by updating the relevant data
structures (Line $5-9$ in Algorithm 3).

  * •

If an interval $I_{t^{\prime\prime}}$ with $l_{t^{\prime\prime}}=t$, then
color for $I_{t^{\prime\prime}}$ is computed by calling Algorithm 1 (Line $13$
in Algorithm 3) and Algorithm 2 (Line $14$ in Algorithm 3).

In the next lemma we prove that $\sf DSusingBPC$ is $3$-competitive similar as
of [9].

1 if  _$I_{t}\neq NULL$_ then

2 $\sf insertMap(\mathcal{I},<t,I_{t}>)$

3 $\sf insertHeap(\sf H^{min}_{L},(l_{t},t))$

4 $\sf insertHeap(\sf H^{min}_{R},(r_{t},t))$

5

6$r_{t^{\prime\prime}}\leftarrow\sf extractMin(\sf H^{min}_{R})$

7 if  _$r_{t^{\prime\prime}}==t$ _ then

8 $\sf insertHeap(\sf H^{min}_{F{\sf idNumber}},{\sf
idNumber}(I_{t^{\prime\prime}}))$

9 $\sf deleteBST(\sf T_{U{\sf idNumber}},{\sf
idNumber}(I_{t^{\prime\prime}}))$

10 $\sf deleteHeap(\sf H^{min}_{R})$

11

12else

13 $l_{t^{\prime\prime}}\leftarrow\sf extractMin(\sf H^{min}_{L})$

14 if  _$l_{t^{\prime\prime}}==t$_ then

15 compute-idNumber$\sf(I_{t^{\prime\prime}})$

16 compute-binNumber-Nextfit$\sf(I_{t^{\prime\prime}})$

17 $\sf insertMap({\sf Color},<t^{\prime\prime},$ $({\sf
idNumber}(l_{t^{\prime\prime}}),$ ${\sf binNumber}(l_{t^{\prime\prime}}))>)$

18 $\sf deleteHeap(\sf H^{min}_{L})$

19

Algorithm 3 $\sf DSusingBPC(I_{t}=[l_{t},r_{t}])$

###### Lemma 6

$\sf DSusingBPC$ is $3$-competitive.

###### Proof

Let $\sf L_{i}$ be the number of bins, each of capacity $B$, created
corresponding to ${\sf idNumber}=i$. Let ${\sf cost}(L_{i})$ correspond to the
summation of the costs of all the intervals whose color consists of ${\sf
idNumber}=i$ and ${\sf binNumber}=L_{i}$. Notice that in the next-fit
strategy, for any two consecutive bins, the sum of the occupied capacity in
the bins is strictly greater than $B$. Therefore, it is easy to observe that
$\frac{{\sf cost}(L_{i})}{B}\geq\frac{L_{i}-1}{2}$ or $L_{i}\leq
2\cdot\frac{{\sf cost}(L_{i})}{B}+1$. Let ${\sf ALG}$ be the total number of
different colors used by $\sf DSusingBPC$. Let $\sf OPT_{off}^{B}$ be the
number of different colors used by an optimal offline algorithm with budget
constraint $B$ for every color. Then,

${\sf ALG}=\sum_{i=1}^{{\sf gIdNumber}}{L_{i}}\leq
2\cdot{\left(\sum_{i=1}^{{\sf gIdNumber}}{\frac{{\sf
cost}(L_{i})}{B}}\right)}+{\sf gIdNumber}\\\ $

It is easy to observe that $\sf OPT_{off}^{B}\geq\sum_{i=1}^{{\sf
gIdNumber}}{\frac{{\sf cost}(L_{i})}{B}}$. In Lemma 3, we proved that ${\sf
gIdNumber}\leq\sf OPT_{off}$, where $\sf OPT_{off}$ is the number of different
colors used by an optimal offline coloring algorithm without any budget
constraints on the colors. Again, it is easy to observe that $\sf
OPT_{off}\leq\sf OPT_{off}^{B}$. Therefore,

${\sf ALG}\leq 2\cdot{\left(\sum_{i=1}^{{\sf gIdNumber}}{\frac{{\sf
cost}(L_{i})}{B}}\right)}+{\sf gIdNumber}\leq 2\cdot\sf OPT_{off}^{B}+\sf
OPT_{off}\leq 3\cdot\sf OPT_{off}^{B}\\\ $

Hence, algorithm $\sf DSusingBPC$ is $3$-competitive.

###### Lemma 7

The update time of $\sf DSusingBPC$ is $O(\log n)$ in the worst-case.

###### Proof

Update time of $\sf DSusingBPC$ is dominated by standard operations in the
heap, binary search tree and map, every operation taking $O(\log n)$ time in
the worst-case. From Lemma 2 and Lemma 5, computing color for an interval
takes $O(\log n)$ time in the worst-case. Therefore, update time $\sf
DSusingBPC$ is $O(\log n)$ in the worst-case.

Using Lemma 6 and Lemma 7, we get the following result.

###### Theorem 2.1

There exists a deterministic $3$-competitive algorithm for online drone
scheduling problem that performs every drone allotment in $O(\log n)$ worst-
case time, where $n$ is the maximum number of requests being served by some
drone at any instance (active requests) through out the handling of the
request sequence.

###  2.4 Improved algorithm using first-fit for computing ${\sf binNumber}$

In this section, we show that the competitive ratio of our result in Theorem
2.1 can be improved by using the first-fit strategy from online bin packing to
compute ${\sf binNumber}$ for an interval.

###### Definition 3

(first-fit) The first-fit algorithm for online bin packing maintains a
collection of active bins in sorted order (increasing order) with respect to
their appearance. When a new item appears, among all the active bins that can
accommodate the new item, it finds the first bin where the item can fit. If no
such bin exists, a new active bin is created, and the item is placed in the
new bin. When the remaining capacity of an active bin becomes zero, the bin is
closed and marked as inactive.

1 if  _$root.max <{\sf cost}(I_{t})$_ then

2 $\sf gBinNumber[{\sf idNumber}(I_{t})]\leftarrow\sf gBinNumber[{\sf
idNumber}(I_{t})]+1$

3 $insertBST(\sf T_{{\sf binNumber}},\sf gBinNumber[{\sf
idNumber}(I_{t})],B-{\sf cost}(I_{t}))$

4 return $\sf gBinNumber[{\sf idNumber}(I_{t})]$

5

6else

7 if  _$root.rem\geq{\sf cost}(I_{t})$_ then

8 if  _$root.left.max <{\sf cost}(I_{t})$ or $root.left=NULL$_ then

9 $decreaseKey(\sf T_{{\sf binNumber}},root,{\sf cost}(I_{t}))$

10 return $root.bin$

11

12 else

13 return compute-binNumber-Firstfit$(root.left,I_{t})$

14

15 else

16 if  _$root.left\neq NULL$ and $root.left.max\geq{\sf cost}(I_{t})$_ then

17 return compute-binNumber-Firstfit$(root.left,I_{t})$

18

19 else

20 return compute-binNumber-Firstfit$(root.right,I_{t})$

21

Algorithm 4 compute-binNumber-Firstfit$(root,I_{t})$

To implement the first-fit strategy to compute ${\sf binNumber}$, for every
${\sf idNumber}$ in the map $\sf idTobinNumber$, we maintain the set of active
bin numbers in a balanced binary search tree $\sf T_{{\sf binNumber}}$ with
the bin number as the key. We store another map $\sf gBinNumber$ with ${\sf
idNumber}$ as key to store the number of bins created corresponding to this
${\sf idNumber}$. Each $node$ in the tree $\sf T_{{\sf binNumber}}$ has five
attributes, (1) $node.bin$: bin number (2) $node.rem$: remaining capacity of
the drone (3) $node.max$: maximum remaining capacity among all the nodes in
the sub-tree rooted at $node$ (4) $node.left$: left child of the $node$ (5)
$node.right$: right child of the node. $insertBST(T,i,c)$ is used to insert a
$node$ in the tree $T$ with the key $i$ and $node.rem=c$.
$decreaseKey(T,node,c)$ is used to decrease the $node.rem$ by $c$ in the tree
$T$. For both operations, we also need to update the $max$ attribute of each
node on the path of $root$ to $node$.

Suppose we need to compute ${\sf binNumber}$ for interval $I_{t}$ with cost
${\sf cost}(I_{t})$ and ${\sf idNumber}(I_{t})$. Let $\sf T_{{\sf binNumber}}$
be the tree stored in the map $\sf idTobinNumber$ corresponding to ${\sf
idNumber}(I_{t})$. If $\sf T_{{\sf binNumber}}$ is empty, then a new tree is
created consisting of a single $node$ with $node.bin=1$,
$node.rem=node.max=B-{\sf cost}(I_{t})$, and both the children set to $NULL$.
$node$ set as the root of this tree and $\sf gBinNumber[{\sf
idNumber}(I_{t})]=1$. Otherwise, Algorithm 4 is called with the $root$ of $\sf
T_{{\sf binNumber}}$ and the interval $I_{t}$ as its’ parameter. The returned
value is set as the bin number of $I_{t}$ i.e., ${\sf binNumber}(I_{t})$. If
$root.max<cost(I_{t})$, then all the nodes have remaining battery capacity
lesser than ${\sf cost}(I_{t})$. For this case, we need to introduce a new bin
by increasing $\sf gBinNumber[{\sf idNumber}(I_{t})]$ by one and set it as the
bin number of $I_{t}$. A $node$ corresponding to this newly introduced bin
needs to be added to the tree with $node.bin=\sf gBinNumber[{\sf
idNumber}(I_{t})]$ and $node.rem=B-{\sf cost}(I_{t})$ (Line $1-4$). Otherwise
( $root.max\geq{\sf cost}(I_{t})$), we need to find a first bin to fit the
interval $I_{t}$. If $root.rem\geq{\sf cost}(I_{t})$ (Line $6$), then either
the root node or a node in its’ left subtree is the first node to accommodate
$I_{t}$. If the maximum remaining capacity of the roots’ left subtree is
lesser than ${\sf cost}(I_{t})$ (Line $7$), then the bin number to the root
node is set as the bin number of $I_{t}$ (Line $9$). Also, we need to decrease
$root.rem$ by ${\sf cost}(I_{t})$ (Line $8$). If $root.left.max\geq{\sf
cost}(I_{t})$, then we need to find the bin into the left subtree of the root
by calling the Algorithm 4 for the $root.left$ node (Line $10-11$ and
$13-14$). If $root.rem<{\sf cost}(I_{t})$ and $root.left.max<{\sf
cost}(I_{t})$, then the bin to fit the interval $I_{t}$ must be in the right
subtree of root. So we call the Algorithm 4 for this case as $root.right$ as
its’ parameter (Line $15-16$). All these operations are standard operations in
a balanced binary search tree and take $O(\log n)$ time in the worst-case.
Therefore, the update time of the $\sf DSusingBPC$ (Algorithm 3) is $O(\log
n)$ in the worst-case if the first-fit strategy is used to compute ${\sf
binNumber}$ in Algorithm 2 instead of next-fit. Here, $n$ is the number of
requests received at any instance.

###### Theorem 2.2

Algorithm $\sf DSusingBPC$ is $2.7$-competitive when the first-fit strategy is
used to compute ${\sf binNumber}$.

###### Proof

Let $\sf L_{i}$ be the number of bins, each of capacity $B$, created
corresponding to ${\sf idNumber}=i$. We apply the same weight function
described in [11] for analyzing the first fit strategy on $L_{i}$. After
applying the weight function, let $\sf W(i)$ correspond to the summation of
the weights of all the intervals whose color consists of ${\sf idNumber}=i$
and ${\sf binNumber}$ is one of the $\sf L_{i}$ bins created corresponding to
${\sf idNumber}=i$. It is shown in [11] that $L_{i}\leq W(i)+1$. Let ${\sf
ALG_{f}}$ be the total number of different colors used by algorithm $\sf
DSusingBPC$ when first-fit is used to compute ${\sf binNumber}$. Let $\sf
OPT_{off}^{B}$ be the number of different colors used by an optimal offline
algorithm with budget constraint $B$ for every color. Thus,

${\sf ALG_{f}}=\sum_{i=1}^{{\sf gIdNumber}}{L_{i}}\leq\sum_{i=1}^{{\sf
gIdNumber}}{\left(W(i)+1\right)}=\sum_{i=1}^{{\sf gIdNumber}}{W(i)}+{\sf
gIdNumber}\\\ $

From the first-fit analysis of [11], $\sum_{i=1}^{{\sf gIdNumber}}{W(i)}\leq
1.7\cdot\sf OPT_{off}^{B}$. In Lemma 3, we proved that ${\sf gIdNumber}\leq\sf
OPT_{off}$, where $\sf OPT_{off}$ is the number of different colors used by an
optimal offline coloring algorithm without any budget constraints on the
colors. Again, it is easy to observe that $\sf OPT_{off}\leq\sf
OPT_{off}^{B}$. Therefore,

${\sf ALG_{f}}\leq{\sf gIdNumber}+\sum_{i=1}^{{\sf gIdNumber}}{W(i)}\leq
1.7\sf OPT_{off}^{B}+\sf OPT_{off}\leq 2.7\sf OPT_{off}^{B}\\\ $

Hence, algorithm $\sf DSusingBPC$ when the first-fit strategy is used to
compute ${\sf binNumber}$ is $2.7$-competitive.

##  3 Online Variable-size Drone Scheduling

In this section, we discuss a variant of the drone scheduling problem, and we
call it online variable-size drone scheduling (OVDS). In the OVDS problem, we
know all the customer requests in advance. So, we compute the underlying
interval graph and partition the intervals into independent sets before the
drones appear online. Each interval is associated with a cost reflecting the
drone’s battery consumption for the corresponding customer request. Drones
with different battery capacities appear online. We allocate as many requests
as possible to the current drone from one of the pre-computed independent
sets. The battery capacity of every drone is at least the maximum cost of an
interval. The goal is to use a minimum number of drones.

Our motivation for the OVDS problem comes from the following scenario: there
is an allocation unit for drones that allocates a drone as and when requested
by a delivery unit. The delivery unit is not aware of the battery capacities
of the drones in advance. Therefore, upon request, when a drone is allotted,
depending upon the drone’s battery capacity, the delivery unit aims to
schedule as many requests as possible for that drone. The delivery unit keeps
raising requests to the allocation unit until all customer requests are
allocated to some drones. The goal is to minimize the number of requests
raised by the delivery unit.

Let $\mathcal{I}$ be the set of intervals corresponding to the customer
requests. To allocate drones to the intervals, we assign a color to the
intervals, similar to Section 2. The color of an interval $I$ denoted by ${\sf
Color}(I)$ corresponds to the drone allotted to serve the corresponding
request. The ${\sf Color}(I)$ is defined as two tuples ${\sf Color}(I)=({\sf
idNumber}(I),{\sf binNumber}(I))$. We can compute the ${\sf idNumber}$ of any
interval $I$ by using Algorithm 1 and the set of intervals with a particular
${\sf idNumber}$ forms an independent set. Let $\\{1,2,\cdots,g\\}$ be the set
of ${\sf idNumber}$. There is a list of intervals corresponding to each of the
${\sf idNumber}$ in $\\{1,2,\cdots,g\\}$ and every list is an independent set.
For each such list, we use the any-fit strategy from the bin-packing problem
and compute the ${\sf binNumber}$ for the intervals in that list. We present
the algorithm $\sf VariableSizeDS$ that computes a color (allocates a drone)
for each interval.

###### Definition 4

(any-fit) A bin packing algorithm is called any-fit algorithm if, after
packing a bin, there is no item left in the list that could still fit into the
bin.

1 for  _each $I\in\mathcal{I}$_ do

2 compute-idNumber$(I)$

3

4$g$ is the total ${\sf idNumber}$ used

5 for  _$i=1$ to $g$_ do

6 Create a list $List_{i}$ and $List_{i}\leftarrow\emptyset$

7

8for  _each $I\in\mathcal{I}$_ do

9 Insert $I$ into $List_{{\sf idNumber}(I)}$

10

11for  _each $i=1$ to $g$_ do

12 Sort the intervals in $List_{i}$ as per the non-decreasing order of their
costs

13

14for  _each $i=1$ to $g$_ do

15 $L_{i}\leftarrow 0$

16 while  _$List_{i}$ is non-empty_ do

17 Request for a new drone

18 $b$ is the capacity of this newly arrived drone

19 $L_{i}\leftarrow L_{i}+1$; $rem\leftarrow b$

20 for  _each $I\in List_{i}$ as per the sorted order_ do

21 if  _${\sf cost}(I)\leq rem$_ then

22 ${\sf binNumber}(I)\leftarrow L_{i}$

23 $rem\leftarrow$ $rem-{\sf cost}(I)$

24 Delete $I$ from $List_{i}$

25

26 else

27 break

Algorithm 5 $\sf VariableSizeDS(\mathcal{I})$

###### Lemma 8

The running time of $\sf VariableSizeDS$ is $O(n\log n)$, where $n$ is the
total number of requests or delivery.

###### Proof

We can compute the ${\sf idNumber}$ of each interval in $O(n\log n)$ time (
follows from Lemma 2). Creating a list for each ${\sf idNumber}$ and assigning
the intervals takes linear time. Thereafter sorting the intervals in each list
takes $O(n\log n)$ time. At last, finding ${\sf binNumber}$ for each interval
is done by some constant comparisons and assignments, which takes linear time
again. Thus, the overall time complexity is $O(n\log n)$ time.

###### Theorem 3.1

$\sf VariableSizeDS$ is $(2\alpha+1)$-competitive, where $\alpha$ is the ratio
between the maximum and minimum battery capacity.

###### Proof

Let $L_{i}$ be the number of different bins or ${\sf binNumber}$ created for
${\sf idNumber}=i(1\leq i\leq g)$. Thus, $\\{(i,1),(i,2),\cdots,(i,L_{i})\\}$
are the set of different colors used for the ${\sf idNumber}=i$. For each
color, there is a request for a new drone. So, there is an associated battery
capacity for each color. Let $B_{i}^{j}$ be the battery capacity of the color
$(i,j)$. Let ${\sf cost}(i,j)$ be the sum of all the costs of the intervals
having color $(i,j)$. Then, ${\sf cost}(i,j)+{\sf cost}(i,j+1)>B_{i}^{j}$, for
$1\leq i\leq g$ and $1\leq j\leq L_{i}-1$. So, ${\sf
cost}(i,1)+2\cdot\sum_{j=2}^{L_{i}-1}{\sf cost}(i,j)+{\sf
cost}(i,L_{i})>\sum_{j=1}^{L_{i}-1}B_{i}^{j}$ , implies
$\sum_{i=1}^{g}\sum_{j=1}^{L_{i}-1}B_{i}^{j}<2\cdot\sum_{i=1}^{g}\sum_{j=1}^{L_{i}}{\sf
cost}(i,j)$.

Let $B_{min}=\underset{1\leq i\leq g;1\leq j\leq L_{i}}{\min}B_{i}^{j}$, and
$B_{max}=\underset{1\leq i\leq g;1\leq j\leq L_{i}}{\max}B_{i}^{j}$. Then,
$\alpha=\frac{B_{max}}{B_{min}}$. Let $OPT$ be the optimum colors used for the
problem in the offline version when we have the knowledge of all the intervals
along with their costs as well as all the drones with different battery
capacities. Then, $OPT\geq\frac{\sum_{i=1}^{g}\sum_{j=1}^{L_{i}}{\sf
cost}(i,j)}{B_{max}}$. Also, $OPT\geq g$, from Lemma 3. Thus, the total number
of colors used by the algorithm is
$\sum_{i=1}^{g}L_{i}=\sum_{i=1}^{g}(L_{i}-1)+g\leq\sum_{i=1}^{g}\sum_{j=1}^{L_{i}-1}(\frac{B_{i}^{j}}{B_{min}})+g$.
Therefore,

$\sum_{i=1}^{g}L_{i}\leq\frac{2\cdot\sum\limits_{i=1}\limits^{g}\sum\limits_{j=1}\limits^{L_{i}}{\sf
cost}(i,j)}{B_{min}}+g\leq\frac{2\cdot\alpha\cdot\sum\limits_{i=1}\limits^{g}\sum\limits_{j=1}\limits^{L_{i}}{\sf
cost}(i,j)}{B_{max}}+g\leq(2\alpha+1)OPT.\\\ $

Hence the proof.

##  4 Interval generator

In Section 2, we studied the online drone Scheduling problem where intervals
appear online with a cost and every interval corresponds to a request of a
customer. Here, we analyze further how these intervals can be generated. For
that reason, we present and study the Interval Generator problem. Let us
consider the following scenario. There is a truck that carries some drones and
moves on a specific route throughout the day. Moreover, on the route have been
predefined stops for the truck where the drones can take off or/and land.
During the movement of the truck, we receive a request from a customer that we
want to serve. The drones will be used to serve the customers. The goal here
is to find the stops that a drone needs to take off and land that minimize a
cost function. Since in the drone scheduling problem for every interval
(request) we want to assign a unique drone, we need to represent the takeoff
and landing time for a particular request as an interval. Thus, the interval
generator problem receives a request in a specific time and produces an
interval for this request with a cost based on the stops and the position of
the truck.

More specifically, a request $q$ is given as coordinates in the plane and
denoted by $q=(x_{q},y_{q})$. Now, for the truck we need to consider two
things, its route and its position on that. The route can be considered as
known and contains some stops. Let
$\mathcal{S}=\\{S_{1},S_{2},\cdots,S_{k}\\}$ be the set of stops on the route.
We assume that the stops $S_{1},\cdots,S_{k}$ are increasingly ordered w.r.t.
the time of visit by the truck. Thus, the stop $S_{1}$ is the start point and
$S_{k}$ is the end point of the truck. Every stop is given as coordinates in
the plane, $S_{i}=(x_{i},y_{i})$ for all $i\in[1,k]$. Additionally, the stops
will be used for the takeoff and landing of the drones as well as for the
loading of the drones with the shipment. Moreover, the position $p$ of the
track corresponds to some coordinates $(x_{p},y_{p})$ on the plane.

Next, the cost of a request can be considered as the time that a drone needs
to serve a customer and to return to the truck (as we assume that the battery
consumption is a linear function of the time). For simplicity, we can assume
that there is no battery consumption during the loading of the drones.
Furthermore, we say that the cost of a request is _invalid_ if the landing
stop of the drone is ‘smaller’ than the position of the truck at the time when
the drone returns. Otherwise, we say that the cost of a request is _valid_. We
define it formally later. Now we can construct an interval for the drone
scheduling problem. In Figure 1, we can see an example of the above scenario
where the route is a straight line. Formally, we have the following definition
for the problem.

###### Definition 5

(Interval Generator)
Input: A set of stops $\mathcal{S}=\\{S_{1},S_{2},\cdots,S_{k}\\}$, where each
$S_{i}=(x_{i},y_{i})$, the truck position $p=(x_{p},y_{p})$, a request
$q=(x_{q},y_{q})$ and a cost function ${\sf cost}(\cdot)$.
Output: An interval $I_{q}=[\ell_{q},r_{q}]$ for the interval graph
$\mathcal{I}$ with the minimum valid cost.

Observe that in the Interval Generator problem there is no restriction on the
speed of the truck and the drones. Thus, we provide some assumptions to
simplify the study of the problem.

Assumptions (1) The speed of the drones and the truck is known and fixed. (2)
The time of takeoff and landing as well as the time of loading and unloading
is considered negligible. (3) The time that the truck remains on a stop is
also considered negligible. (4) The timeline of the day is known.

Figure 1: Above: The route of the truck is a straight line. The stops are
given as black solid dots, the position $p$ of the truck is given as a blue
square on the route and the request $q$ as a red square. Moreover, $S_{i}$ is
the takeoff stop and $S_{j}$ is the landing stop and $I_{q}=[\ell_{q},r_{q}]$
is the corresponding interval of the solution. Below: The timeline of the day
with the corresponding times $T_{1},\cdots,T_{k}$ of the stops. $T_{q}$ is the
position of the truck in the timeline when we receive the request $q$.

We describe now the way that we create an interval for a request for the
Interval Generator problem. This can be done in two steps. First, we find the
stops that minimize the cost of the request and after we compute the
corresponding interval. For the first part, by the input of the problem we
have a set of stops $S$ and based on the assumptions we know the location of
the truck every moment through the day. Hence, we can correspond every stop
$S_{i}$ to a specific time $T_{i}$ in the timeline of the day. Moreover, when
we receive a request $q$ we know the position $p$ of the truck and thus we can
correspond the position of the truck as a time $T_{q}$ in the timeline in
order to know when we have received the request. Next, we need to find the
stops (times) for the takeoff and landing that minimize the cost of the
request.

As mentioned earlier the cost of a request is the total time from the takeoff
stop to the position of the request and from there to the landing stop. Assume
$S_{i}$ is the takeoff stop and $S_{j}$ is the landing stop; then we denote
the total time of the delivery as $(S_{i},q,S_{j})$. Thus, in order to find
the minimum cost we need to try all possible pairs of stops with $i<j$. Also,
the cost of the request needs to be valid. This means that the time of
$(S_{i},q,S_{j})$ is not greater than the time that the truck needs to go from
the stop $S_{i}$ to the stop $S_{j}$. Let $A$ be the set of all pairs that
give valid cost. Hence, the formula that computes the minimum valid cost for a
request $q$ is the following: ${\sf cost}(q)=\min_{p\leq i<j,\,(i,j)\in
A}(S_{i},q,S_{j}).$

So far we have found the optimal stops for the delivery. Let $S_{i}$ and
$S_{j}$ be the takeoff and landing stop respectively. We have corresponded
every stop to a specific time in the timeline of the day. Thus, we know the
takeoff time $T_{i}$ and the landing time $T_{j}$. Now we are ready to create
an interval $I_{q}=[\ell_{q},r_{q}]$ with left endpoint $\ell_{q}$, right
endpoint $r_{q}$ and to assign a cost to the interval. By the definition of
the interval graphs, every vertex corresponds to an interval of the real line
in such a way that two vertices are adjacent if and only if the corresponding
intervals have a nonempty intersection. We can consider the timeline of the
day as the real line and the times $T_{i}$ as points on the real line. Then
every valid delivery corresponds to an interval for the interval graph with
left endpoint $\ell_{q}=T_{i}$ and right endpoint $r_{q}=T_{j}$. Moreover, we
associate the cost of the request $q$ to this interval
$I_{q}=[\ell_{q},r_{q}]$, that is ${\sf cost}(I_{q})={\sf cost}(q)$. Hence, we
have described how we can find a solution for the Interval Generator problem.
The Algorithm 6 computes an optimal solution for the Interval Generator
problem.

1 Input: Stops $\mathcal{S}=\\{S_{1},S_{2},\cdots,S_{k}\\}$, truck position
$p=(x_{p},y_{p})$, request $q=(x_{q},y_{q})$ and cost function ${\sf
cost}(\cdot)$

2 Output: An interval $I_{q}=[\ell_{q},r_{q}]$ with the minimum valid cost.

3 Let $M$ be the set of the indexes of the stops that are $\geq p$.

4 for  _ $i\in M$_ do

5 for  _ $j\in M\quad and\quad i<j$_ do

6 if  _ $(i,j)$ is a valid pair_ then

7 if  _${\sf cost}(q) >(S_{i},q,S_{j})$_ then

8 ${\sf cost}(q)=(S_{i},q,S_{j})$ and $(i^{\prime},j^{\prime})=(i,j)$

9

10Find the times $T_{i^{\prime}}$ and $T_{j^{\prime}}$

11 Set $\ell_{q}=T_{i^{\prime}}$ and $r_{q}=T_{j^{\prime}}$

12 Set ${\sf cost}(I_{q})={\sf cost}(q)$

13 Return $I_{q}=[\ell_{q},r_{q}]$

14

Algorithm 6 $\sf IntervalGenerator$

###### Theorem 4.1

The running time of the Algorithm 6 is $O(k^{2}\cdot\mathcal{T})$, where
$\mathcal{T}$ is the time needed to compute each of $(S_{i},q,S_{j})$.

###### Proof

We know that the set $M$ contains at most $k$ elements. Since we have two
loops over this set, we have at most $k^{2}$ iterations. In each iteration, we
need first to compute the cost $(S_{i},q,S_{j})$. This can be done in
$\mathcal{T}$ time. Next, we check if the cost is valid and if it is we
compare it with the previous minimum cost. This can be done in constant time
since the time that the truck needs to move from the stop $S_{i}$ to the stop
$S_{j}$ is known by the assumptions. Thus, the running time of the Algorithm 6
is $O(k^{2}\cdot\mathcal{T})$.

In Section 2, we assume that all intervals have distinct endpoints or boundary
points. However, Algorithm 6 can generate multiple intervals with common
boundary points. To make each boundary point distinct, a small $\epsilon>0$
can be used for such intervals. For example $I_{q}=[\ell_{q},r_{q}]$ and
$I_{p}=[\ell_{p},r_{p}]$ with one common boundary point, i.e.,
$r_{q}=\ell_{p}$. These can be transformed into $[\ell_{q},r_{q}-\epsilon]$
and $[\ell_{p}+\epsilon,r_{p}]$ to ensure all boundary points of the $I_{p}$
and $I_{q}$ are distinct.

##  5 Conclusion

In this paper, we have studied a hybrid truck-drone model for last-mile
delivery where delivery requests appear online during the moving truck
carrying multiple drones in a predefined path. A drone is enabled to launch
from the truck to deliver a parcel to the customer, as per the request, and
after delivery, it returns and lands on the truck. The objective is to
minimize the number of drones allocated to complete all feasible deliveries.
To address this, we have proposed two algorithms for online drone scheduling
problem based on next-fit and first-fit bin-packing strategies. We have
analyzed the performances of the algorithms in terms of competitive ratio and
worst-case time complexity. We also have introduced online variable-size drone
scheduling problem and proposed an algorithm to solve it. We have discussed
the interval generator, a tool to calculate delivery time intervals for each
request. The interval generator aims to minimize the delivery cost associated
with each request, choosing delivery intervals from predefined discrete
stopping points on the truck’s path. We can treat these points as mini-
warehouses where packages can be loaded into the truck, if necessary, to
fulfill the requests. For future work, evaluating the optimal drone-truck path
to deliver all types of packages would be interesting along with a simulation
study for possible practical application. Another interesting problem to
address in the future is determining the delivery interval for each request
when there are no predefined discrete stopping points along the truck’s path.
The model’s applicability can also be extended to a drone-delivery system with
a recharging battery policy for the drones.

## References

  * [1] Ahani, G., Yuan, D., Zhao, Y.: Age-optimal uav scheduling for data collection with battery recharging. IEEE Communications Letters 25(4), 1254–1258 (2020)
  * [2] Amazon: Amazon customers in lockeford, california, will be among the first to receive prime air drone deliveries in the u.s. (2022), https://www.aboutamazon.com/news/transportation/amazon-prime-air-prepares-for-drone-deliveries
  * [3] Betti Sorbelli, F., Corò, F., Das, S.K., Palazzetti, L., Pinotti, C.M.: On the scheduling of conflictual deliveries in a last-mile delivery scenario with truck-carried drones. Pervasive and Mobile Computing 87, 101700 (2022)
  * [4] Boysen, N., Briskorn, D., Fedtke, S., Schwerdfeger, S.: Drone delivery from trucks: Drone scheduling for given truck routes. Networks 72, 506 – 527 (2018)
  * [5] Chiang, W.C., Li, Y., Shang, J., Urban, T.L.: Impact of drone delivery on sustainability and cost: Realizing the uav potential through vehicle routing optimization. Applied Energy 242, 1164–1175 (2019)
  * [6] Crişan, G.C., Nechita, E.: On a cooperative truck-and-drone delivery system. Procedia Computer Science 159, 38–47 (2019), knowledge-Based and Intelligent Information & Engineering Systems: Proceedings of the 23rd International Conference KES2019
  * [7] Daknama, R., Kraus, E.: Vehicle routing with drones. ArXiv abs/1705.06431 (2017)
  * [8] DHL: Dhl express launches its first regular fully-automated and intelligent urban drone delivery service (2019), https://www.dhl.com/global-en/home/press/press-archive/2019/dhl-express-launches-its-first-regular-fully-automated-and/intelligenturban-drone-delivery-service.html
  * [9] Epstein, L., Levin, A.: On bin packing with conflicts. SIAM J. Optim. 19(3), 1270–1298 (2008), https://doi.org/10.1137/060666329
  * [10] Fulkerson, D.R., Gross, O.A.: Incidence matrices and interval graphs. Pacific Journal of Mathematics 15, 835–855 (1965)
  * [11] Garey, M.R., Graham, R.L., Johnson, D.S., Yao, A.C.: Resource constrained scheduling as generalized bin packing. J. Comb. Theory, Ser. A 21(3), 257–298 (1976), https://doi.org/10.1016/0097-3165(76)90001-7
  * [12] Ghazzai, H., Menouar, H., Kadri, A., Massoud, Y.: Future uav-based its: A comprehensive scheduling framework. IEEE Access 7, 75678–75695 (2019)
  * [13] Goodchild, A., Toy, J.: Delivery by drone: An evaluation of unmanned aerial vehicle technology in reducing co2 emissions in the delivery service industry. Transportation Research Part D: Transport and Environment 61, 58–67 (2018)
  * [14] Jaillet, P., Wagner, M.R.: Online vehicle routing problems: A survey. The Vehicle Routing Problem: Latest Advances and New Challenges pp. 221–237 (2008)
  * [15] Jana, S., Mandal, P.S.: Approximation algorithms for drone delivery packing problem. In: Proceedings of the 24th International Conference on Distributed Computing and Networking. pp. 262–269 (2023)
  * [16] Jana, S., Mandal, P.S.: Approximation algorithms for drone delivery scheduling with a fixed number of drones. Theoretical Computer Science 991, 114442 (2024)
  * [17] Kim, S.J., Lim, G.J.: A hybrid battery charging approach for drone-aided border surveillance scheduling. Drones 2(4), 38 (2018)
  * [18] Kim, S., Kim, T., Suh, K., Jeon, J.: Energy and environmental performance of a passenger drone for an urban air mobility (uam) policy with 3d spatial information in seoul. Journal of Cleaner Production 415, 137683 (2023)
  * [19] Laporte, G.: The traveling salesman problem: An overview of exact and approximate algorithms. European Journal of Operational Research 59(2), 231–247 (1992)
  * [20] Laporte, G.: The vehicle routing problem: An overview of exact and approximate algorithms. European journal of operational research 59(3), 345–358 (1992)
  * [21] Murray, C.C., Chu, A.G.: The flying sidekick traveling salesman problem: Optimization of drone-assisted parcel delivery. Transportation Research Part C: Emerging Technologies 54, 86–109 (2015)
  * [22] Murray, C.C., Raj, R.: The multiple flying sidekicks traveling salesman problem: Parcel delivery with multiple drones. Transportation Research Part C: Emerging Technologies 110, 368–398 (2020)
  * [23] Park, S., Zhang, L., Chakraborty, S.: Battery assignment and scheduling for drone delivery businesses. In: 2017 IEEE/ACM International Symposium on Low Power Electronics and Design (ISLPED). pp. 1–6. IEEE (2017)
  * [24] Wang, X., Poikonen, S., Golden, B.: The vehicle routing problem with drones: several worst-case results. Optimization Letters 11, 679–697 (2017)