#

Francesco Betti Sorbelli  Department of Computer Science and Mathematics,
University of Perugia, Italy  Federico Corò  Department of Computer Science
and Mathematics, University of Perugia, Italy  Sajal K. Das  Department of CS,
Missouri University of Science and Technology, USA  Cristina M. Pinotti
Department of Computer Science and Mathematics, University of Perugia, Italy
Anil Shende  Department of Math., CS, & Physics, Roanoke College, USA

###### Abstract

In this paper, we present a drone-based delivery system that assumes to deal
with two different mixed-areas, i.e., rural and urban. In these mixed-areas,
called EM-grids, the distances are measured with two different metrics, and
the shortest path between two destinations concatenates the Euclidean and
Manhattan metrics. Due to payload constraints, the drone serves a single
customer at a time returning back to the dispatching point (DP) after each
delivery to load a new parcel for the next customer. In this paper, we present
the $1$-Median Euclidean-Manhattan grid Problem (MEMP) for EM-grids, whose
goal is to determine the drone’s DP position that minimizes the sum of the
distances between all the locations to be served and the point itself. We
study the MEMP on two different scenarios, i.e., one in which all the
customers in the area need to be served (full-grid) and another one where only
a subset of these must be served (partial-grid). For the full-grid scenario we
devise optimal, approximation, and heuristic algorithms, while for the
partial-grid scenario we devise optimal and heuristic algorithms. Eventually,
we comprehensively evaluate our algorithms on generated synthetic and quasi-
real data.

##  1 Introduction

Drones or Unmanned Aerial Vehicles (UAVs) are recently becoming widely used in
civil applications such as environmental protection [1, 2, 3, 4], public
safety [5, 6, 7], localization [8], and smart agriculture [9, 10]. Currently,
there is a growing interest in the use of drones in smart cities [11]. This
interest is particularly increased after the global presence of the COVID-19
disease [12, 13]. Recently, [14] discuss about the use of drones as carriers
for distributing and transporting drugs and medicines, and revealed that 86.7%
of people agree that drones are more effective, surely faster, and less
polluting than any ground-based distribution system. Quarantine, closure of
borders, and social distancing forced people to stay indoors for long periods,
allowing them to go out only for essential activities. Therefore, considering
the need to avoid all unnecessary direct human contact, people started to rely
heavily on online stores for their regular shopping. In parallel, large
companies like Amazon are testing drone-based delivery systems, particularly
for what is known as “last-mile” small item logistics. For example, Amazon
introduced “Amazon Air Prime”, a service that uses drones able to deliver
goods up to $25\mathrm{\,kg}$ to customers within a radius of
$16\mathrm{\,km}$ [15, 16, 17], and Domino’s developed a pizza-delivery
service using drones [18].

There are countless advantages to using drones for deliveries, including
economic benefits, saving on greenhouse gas, and the ability to deliver in
time-critical situations or hard-to-reach places. With the growth of
commercial interest, researchers have begun to study variants of the Traveling
Salesman Problem (TSP) for drones [19]. However, due to the payload constraint
that forces the drone to return after each delivery, TSP is not suitable for
drones [20, 21]. Of particular interest is the work proposed by [22], where a
drone combined with a truck is used to make multiple deliveries in a given
area, going back and forth from the truck. In that case, it is crucial to find
the best location for the truck to minimize the distance traveled by the
drone.

(a) Rural area: Miami, Florida, US.

(b) Urban area: New York City, New York, US.

(c) Mixed area: Chicago, Illinois, US.

Figure 1: Examples of areas. In the rural area the drone freely moves on
straight lines according to the Euclidean metric (a), while in the urban area,
due to the presence of skyscrapers, the drone can only move over the streets
according to the Manhattan metric (b). In the mixed area the drone has to
combine both the metrics (c).

In this article, we imagine offering a drone delivery service to the customers
of a delivery area that covers two mixed areas, i.e., a rural and/or an urban
one. In urban areas (see Figure 1b), for safety and privacy reasons (see
regulation of inhabited centers, e.g., [23], it is assumed that the drone
flies over the streets since it cannot fly beyond a certain maximum allowed
altitude. Namely, in urban areas, it might be forbidden for the drone to
travel along the straight line connecting the dispatching point (DP) and the
delivery if such a straight line, due to tall obstacles, exists only at an
altitude higher than the maximum altitude allowed by the regulations. However,
the drone can at least save time – if not distance traveled – in the urban
area. In fact, a drone is usually faster than a conventional wheeled vehicle
in a crowded area: it will certainly not get stuck in traffic. In rural areas
(see Figure 1a), instead, due to the fact that tall buildings are not present,
the drone can move freely and follow the shortest path to reach its
destination, making clear its advantage on the traveled distance in this case.
In the case of mixed areas (see Figure 1c), i.e., rural and urban, the drone
must follow both the Euclidean and Manhattan metrics. In particular, in the
real example provided in Figure 1c, the Euclidean part resides on the left,
while the Manhattan part resides on the right. Actually, the two areas are
split by the river.

We model the proposed drone delivery area as a two-dimensional mixed-grid. The
vertices of the grid represent both the possible locations of the drone’s DP
and the possible delivery destinations, and they are placed in rows and
columns as in a regular 2-D grid. To model the two different types of areas
through which the drone can move, we have divided the grid into two parts: the
rural area and the urban area. In the former, the distance between two
destinations respects the Euclidean metric, while in the latter, the metric is
Manhattan (taxicab geometry). We are then interested in finding where to set
the DP in order to minimize the sum of the distances between a subset of
delivery destinations and the DP itself. Due to strict payload constraints,
the drone must return back to the DP after each delivery. Therefore, the drone
needs to travel back and forth from the DP as many times as the number of
delivery destinations on the grid. In this paper, we consider two delivery
scenarios. The full-grid scenario where each point of the grid has to be
served by the drone. In light of the COVID-19 pandemic, this scenario happens,
for example, if the drone is used to deliver meals in a lockdown area, or to
deliver self-tests to sick people. In this particular example, the full-grid
scenario is also justified by the fact that after every delivery the drone
should be properly sanitized and disinfected before performing the next
delivery, as proposed by [24]. The partial-grid scenario, instead, assumes
that only a subset of points of the grid is a delivery site. This is the usual
scenario in a delivery system.

The DP selection problem in logistics, while sharing similarities with the
classical facility location problem, requires a more complex modeling approach
that considers factors such as varying demand and transportation costs, and it
involves optimizing multiple metrics simultaneously (Euclidean and Manhattan),
which is not the case for the original facility location problem.

In this paper, we present extensions and improvements to our earlier work
discussed in two conference papers, i.e., [25] and [26]. In [25] we focused on
the full-grid scenario. In this paper we extend our work to include a partial-
grid scenario. We suitably adapt results and ideas from [26] to this scenario,
and provide new efficient algorithms for both the scenarios. The contributions
of this paper are summarized as follows.

  * •

We introduce the EM-grid model, which characterizes the delivery area for the
drone, formed by two contiguous areas, i.e., one rural and one urban, that
follow the Euclidean and Manhattan metrics, respectively.

  * •

We define the $1$-Median Euclidean-Manhattan grid Problem (MEMP) and devise
time-efficient algorithms for the full-grid and partial-grid scenarios. For
the full-grid scenario we devise optimal, approximation, and heuristic
algorithms, while for the partial-grid scenario we devise optimal and
heuristic algorithms. We also give all the proofs for the full-grid scenario
that had been left blank in the previous conference paper [25].

  * •

In addition to comparing the performance of our presented algorithms on
randomly generated synthetic data, we also evaluate their effectiveness on
quasi-real data obtained by adapting real city maps with our proposed grid
model, providing a more realistic assessment of their practical utility.
Furthermore, to ensure the accuracy of our comparison, we incorporated data
from a real drone to measure the distances traveled by the dispatching point
in our evaluation on quasi-real data.

The rest of the paper is organized as follows. Section 2 reviews the related
work. Section 3 formally defines MEMP. Section 4 and Section 5 describe
properties and algorithms for efficiently solving MEMP with full-grid and
partial-grid scenarios, respectively. Section 6 evaluates our algorithms, and
Section 7 offers conclusions and future research directions.

##  2 Related Work

In the literature, many works attempt to solve the drone-based last-mile
delivery problem. To the best of our knowledge, drones have been considered in
a delivery system for the first time by [27]. Specifically, they study the
cooperation between a truck and a drone to deliver packages to customers. The
problem to solve is the Flying Sidekick Traveling Salesman Problem (FSTSP),
which is a variant of the classic TSP. In the FSTSP, a drone can autonomously
perform deliveries to the customers directly flying from the main depot or can
be helped by a truck. In the latter case, the drone flies from the truck,
delivers the package, and then rendezvouses with the truck again in a third
location. However, when the drone flies, the truck can do other deliveries
independently, but still, it has to wait for the drone at the rendezvous
location. For solving FSTSP, the authors propose an optimal mixed-integer
linear programming (MILP) and two heuristics for solving instances of
practical sizes. Then, [28] investigate the same scenario with multiple drones
in by introducing the Multiple Flying Sidekicks Traveling Salesman Problem
(mFSTSP). Even for the mFSTSP, they provide an optimal MILP formulation along
with a heuristic solution approach that consists of solving a sequence of
three sub-problems.

Recently, [29] present an exact formulation for FSTSP while also simplifying
the model reducing the number of constraints and thus be able to solve several
benchmark instances from the literature. However, in these works, the drones
fly according to the Euclidean metric. [30] introduce the multiple Traveling
Salesman Problem with Drone Stations (mTSP-DS), which is an extension to the
classical multiple Traveling Salesman Problem (mTSP). In this problem,
multiple trucks starting to/from a single depot are in charge of supplying
some packet stations that host autonomous vehicles (drones or robots). On
these stations, each truck can launch and operate drones/robots to serve
customers. The objective of mTSP-DS is to serve all customers minimizing the
makespan. The problem is formulated as an MILP only suitable to solve small
instances. For larger instances, many matheuristic algorithms are presented.
[31] introduce the Vehicle Routing Problem with Drones and En Route Operations
(VRPDERO), which is an extension to the Vehicle Routing Problem with Drones
(VRPD). In this problem, drones may not only be launched and retrieved at
vertices but also on some discrete points that are located on each arc. The
problem is formulated as an MILP, and matheuristic approaches are presented to
deal with large instances. The goal of both [31, 30] is to minimize the
makespan, while ours is to find the best DPs from where to launch the drones.

[25] introduce the drone-based delivery area modeled as EM-grids where a drone
is used for delivering small packages to customers. Given the delivery area
divided into two contiguous areas, i.e., the rural and the urban areas, the
goal is to find the optimal DP (depot or warehouse) for the drone in order to
minimize the sum of all the distances between all the potential customers and
the DP itself. However, due to strict payload constraints, the drone cannot
serve more than a customer at a time, and after each delivery, the drone must
go back to the depot. A similar approach, but in a different context, has been
studied by [26]. In such a scenario, the objective is to determine the optimal
cart point for the drone that minimizes the distances between a set of items
(on shelves) and the cart itself, assuming that shelves follow the two
aforementioned metrics. Differently from [25], in [26] only a subset of points
needs to be considered. Moreover, [26] compare the current human-based system
with respect to the newly proposed one based on drones.

[32] compare and contrast the performance of many algorithms in a delivery
scenario. They model the delivery area as a circular region with a central
depot, while customers are randomly distributed throughout the region.
Different temporal and spatial metrics are compared when evaluating these
algorithms. In particular, they evaluate the impact of having distances
measured according to both the Euclidean and Manhattan distance metrics. The
number of customers stochastically varies under both Manhattan and Euclidean
distance metrics. The paper states that the number of customer deliveries and
the metric used to measure travel distance impacts a decision maker’s choice
of the best algorithm and that employing multiple algorithms is recommended.

[33] explore the implications and advantages of strategic planning on urban
delivery services. More specifically, the preferred method and local impacts
of vehicle trips may vary by neighborhood characteristics (e.g., traffic or
customer demands). Instead of searching for an optimal route, the paper
focuses on the estimation of the vehicles’ miles traveled (VMT) per meal
order, considering different types of neighborhoods, delivery scenarios, and
strategies. The proposed system is tested and evaluated in Chicago, showing
that alternative delivery strategies can greatly reduce the VMT per order
based on the type of neighborhood. In the evaluation, both the Euclidean and
Manhattan metrics are combined. However, although different metrics have been
evaluated, drones have not been used in either [32] or [33].

Recently, there has been more effort in solving the last-mile delivery problem
using drones instead of a standard vehicle, due to the flexibility of drones
and their capability to fly over obstacles and avoid traffic. [34] investigate
the problem of solving the TSP with a Drone (TSP-D) where the drone rides on
the truck. [35] investigate the cooperation between a truck and multiple
drones. Each delivery is characterized by a drone’s energy cost, a reward
based on its priority, and a time interval (launch and rendezvous with the
truck). This work aims at finding an optimal scheduling for the drones that
maximizes the overall reward, subject to the drone’s battery capacity while
ensuring that the same drone performs deliveries that do not overlap. Results
show that the presented problem is $NP$-hard, therefore, different heuristics
for solving the problem in a time-efficient way are proposed. More recently,
[36] investigate the feasibility of performing deliveries with a drone in the
presence of external factors such as wind.

[37] compare the traditional truck-based delivery system against the drone-
based one to reduce the general energy consumption and hence reduce the gas
emissions. They also take into account traffic congestion. They propose a
mixed-integer green routing model with traffic restrictions and a genetic
algorithm to efficiently solve the complex routing problem, showing that
drones can accomplish more deliveries and at a lower cost (in terms of CO2
emissions and energy consumption) compared to standard methods and that
traffic types impact the results.

[38] study the last-mile scenario formed by multiple drones assisted by a
single truck that carries them. The customers to be served by the drones form
a clustering, and each drone is assigned to a specific cluster. The initial
position of a drone is called a “cluster focal point”. Once all the focal
points are computed, the truck needs to visit these points by minimizing its
traveled route. Moreover, due to payload constraints, the drones serve their
customers one at a time. The truck cannot follow the Euclidean metric, while
the drones do. The authors propose an optimal mixed integer nonlinear
programming (MINLP) solution as well as an unsupervised machine learning-based
heuristic algorithm.

[39] present a variation on the theme. The proposed model assumes that drones
are assisted by trucks, that carry them through the city. The trucks start
from the main depot and park the drones in specific locations where drones
have to serve the customers by a sequence of back-and-forth flights. Both the
trucks and the drones move according to the Euclidean metric. The introduced
problem is called Energy Minimizing and Range Constrained Drone Delivery
Problem (ERDDP) whose objective is to minimize the total operational cost
including an explicit calculation of the energy consumption of the drone as a
function of the drone’s speed. The ERDDP is formulated as a second order cone
program instance.

[40] propose a last-mile drone delivery scenario where multiple drones can
exploit charging stations to replenish their batteries. In this setting, the
drones can fly from the main hub to the terminal station and serve, one at a
time, a subset of customers in the neighborhood, and then go back to the hub.
Alternatively, they can move from a terminal station to a charging station to
refill the battery and perform subsequent deliveries to other neighborhoods
belonging to other terminal stations. The objective function aims to either
minimize the number of charging stations or to minimize the overall traveled
distance. The authors solve this problem by proposing an optimal and a
heuristic solution.

[41] combine the pickup and the delivery requests in a system with stations
(nodes in a given graph) and introduce the Hybrid Vehicle-Drone Routing
Problem. Vehicles visit stations to transport delivery items and drones, while
drones are launched and collected only at stations. The problem is formulated
as a mixed-integer program, which minimizes the vehicle and drone routing cost
to serve all customers. To solve the problem, the authors use an extension of
the classic Clarke and Wright algorithm (see [42], a known heuristic to solve
the Vehicle Routing Problem. We remark that while their goal is to find the
best route for drones and trucks, our goal is to find the best stations, i.e.,
the DPs, from where to launch the drones.

Finally, [43] investigate the problem of placing drone charging facilities in
an area to help increase the coverage range of drones for commercial
deliveries. The authors present an MILP formulation, and then a heuristic is
proposed to effectively solve the problem.

##  3 Problem Definition

In this section, we first introduce the delivery area model and how the drone
moves inside it, and then formally describe the problem to solve.

###  3.1 Delivery Area Model

To model the delivery area, we define the Euclidean-Manhattan-Grid (EM-grid)
as $G=(R,C,K)$, that is a 2-D grid with $R$ rows, $C$ columns, and the Border
$B$ is the column $K\in[1,C]$ that separates the Euclidean grid $E$ (rural
area) from the Manhattan grid $M$ (urban area) (see Figure 2). Specifically,
$E=\\{1,\dots,R\\}\times\\{1,\dots,K\\}$,
$B=\\{1,\ldots,R\\}\times\\{K\\}\subseteq E$, and
$M=\\{1,\dots,R\\}\times\\{K+1,\dots,C\\}$.

Figure 2: The EM-grid: $E$ models the rural area (see Figure 1a), while $M$
models the urban area (see Figure 1b).

We assume that the drone delivery system covers a rectangular area: $E$ and
$M$ have the same number of rows $R$. The border consists of a single column,
i.e., $K$. Conventionally, the area consists only of a rural area region if
$K=C$ (i.e., $M=\emptyset$), whereas it is effectively limited to an urban
area if $K=1$. In an EM-grid, there are vertices and edges connecting adjacent
vertices. Any internal vertex $u=(r_{u},c_{u})$ of $G$, i.e., with $1<r_{u}<R$
and $1<c_{u}<C$, is connected to the four adjacent vertices $(r_{u},c_{u}\pm
1)$ and $(r_{u}\pm 1,c_{u})$; whereas, in general, any vertex of the grid,
i.e., with $1\leq r_{u}\leq R$ and $1\leq c_{u}\leq C$, is connected only to
the existing adjacent vertices (i.e., an external vertex has only three or two
neighbors). For simplicity, we assume that the distance between any two pairs
of consecutive vertices on the same row or column is constant and unitary, and
so the weight of any edge is unitary. In this work, the “distance” is a
measure of the time required for performing the delivery, or of the needed
energy for shipping the package. We also assume that every customer can be
reached by the drone to/from the DP. At the DP, the drone can recharge its
battery, or just replace it with fresh ones. Let
$\overline{R}=\lceil\frac{R}{2}\rceil$ be the middle row, let
$\overline{C}=\lceil\frac{C}{2}\rceil$ be the middle column, and let
$\overline{K}=\lceil\frac{K}{2}\rceil$ be the column that halves the Euclidean
sub-grid.

For any two vertices $u$, $v$ in $G$, the distance $d(u,v)$ is the length of
the shortest path traversed by the drone in the EM-grid to go from vertex $u$
to the destination $v$. The Euclidean and Manhattan distances are defined,
respectively, as $d_{E}(u,v)=\sqrt{(r_{u}-r_{v})^{2}+(c_{u}-c_{v})^{2}}$ and
$d_{M}(u,v)=|r_{u}-r_{v}|+|c_{u}-c_{v}|$.

We note that the shortest path between a vertex $u\in E$ and a vertex $v\in M$
is given by $\min_{w\in B}\\{d_{E}(u,w)+d_{M}(w,v)\\}$. In Lemma 1 we prove
that such path is unique and passes through the vertex on the border $B$ that
has the same row as $v$.

######  Lemma 1.

Consider an EM-Grid $G=(R,C,K)$. Given $u=(r_{u},c_{u})\in E$ and
$v=(r_{v},c_{v})\in M$, then $d(u,v)=d_{E}(u,h)+d_{M}(h,v)$ with
$h=(r_{v},K)$.

###### Proof.

Consider the vertex $h=(r_{v},K)$ which shares the same row as $v$, and
another vertex $w=(i,K)$ in $B$ with $h\not=w$. We want to prove that:
$d_{E}(u,h)+d_{M}(h,v)\leq d_{E}(u,w)+d_{M}(w,v)$. This follows by the
triangle inequality applied to the vertices $u,w,h$, i.e., $d_{E}(u,h)\leq
d_{E}(u,w)+|r_{v}-i|=d_{E}(u,w)+d_{M}(w,h)$. ∎

Thus, from now on, $d(u,v)$ is given by:

$\begin{array}[]{ll}d(u,v)=\begin{cases}d_{E}(u,v)&\text{if $u,v\in E$}\\\
d_{M}(u,v)&\text{if $u,v\in(M\cup B)$}\\\ d_{E}(u,h)+d_{M}(h,v)&\text{if $u\in
E,v\in M$ where $h=(r_{v},K)\in B$}\\\ d_{M}(u,h)+d_{E}(h,v)&\text{if $u\in
M,v\in E$ where $h=(r_{u},K)\in B$}\end{cases}\end{array}$ (1)

###  3.2 The Column-Cost

Let the column-cost be the distance traversed by the drone starting from a
vertex $u$ in row $\overline{R}$ to serve all the vertices of a given column,
where the column is on the same side (Euclidean or Manhattan) of the grid as
the vertex $u$. Such column-cost depends on which side of the grid the column
and the vertex $u$ reside and on the number of rows $R$. Let $\Delta_{E}({j})$
be the column-cost to serve a Euclidean column at distance $j$ from the
candidate DP $u=(\overline{R},c_{u})$ with $u\in E$. Similarly, let
$\Delta_{M}({j})$ be the column-cost to serve a Manhattan column at distance
$j$ from the drone’s DP $u=\left(\overline{R},c_{u}\right)$ with $u\in M\cup
B$. One can easily find that:

$\displaystyle\Delta_{E}({j})$
$\displaystyle=j+2\sum\limits_{i=1}^{\overline{R}-1}\sqrt{i^{2}+j^{2}}+((R-1)\bmod
2)\sqrt{\overline{R}^{2}+j^{2}}$ (2) $\displaystyle\Delta_{M}({j})$
$\displaystyle=j+2\sum\limits_{i=1}^{\overline{R}-1}(i+j)+((R-1)\bmod
2)\left(\overline{R}+j\right)$
$\displaystyle=\overline{R}(\overline{R}+1+2j)+j+((R-1)\bmod
2)(\overline{R}+j)$ (3)

###  3.3 The Problem Formulation

In our scenario, the fundamental task is to serve, with the aid of a drone,
the customers of an area, e.g., to distribute viral tests to potentially
infected patients. Due to payload constraints (and, e.g., to avoid the spread
of the disease), the drone cannot serve all the customers on the same flight,
and it has necessarily to go back and forth from a specific position inside
the delivery area (called DP, where the drone, e.g., can be also sanitized
each time) to all the customers. Specifically, this DP is the point in which
all the products for customers are initially stored. Hence, the objective is
to minimize the distance traveled by the drone when it moves inside the
delivery area. We denote this problem as the $1$-Median Euclidean-Manhattan
grid Problem (MEMP) since the goal is to find a single DP inside EM-grids.

Given an EM-grid $G=(R,C,K)$ and a subset of vertices $H\subseteq G$, for an
arbitrary vertex $u=(r_{u},c_{u})\in G$, we define the cost of delivery from
$u$ to each point in $H$, denoted by $\mathcal{C}(H,u)$, as:

$\mathcal{C}(H,u)=2\sum_{v\in H}d(v,u)$ (4)

As noted above, the distance between points $u$ and $v$ is a measure of the
cost of delivery from $u$ to $v$, and the multiplicative constant 2 is in
consideration of the round trip for each delivery. Given $H\subseteq G$, let
$H_{E}$ and $H_{M}$ be the subset of points that lie in the Euclidean and
Manhattan grid, respectively, such that $H_{E}\cup H_{M}=H$ and $H_{E}\cap
H_{M}=\varnothing$. The set $H$ if formed by $n=|H|$ vertices, where
$n_{E}=|H_{E}|$ and $n_{M}=|H_{M}|$, such that $n=n_{E}+n_{M}$.

When $H=G$, Eq. (4) can be rewritten as:

$\mathcal{C}(G,u)=2\sum_{v\in
G}d(v,u)=2\sum_{r_{v}=1}^{R}\sum_{c_{v}=1}^{C}d((r_{v},c_{v}),u)$ (5)

with $v=(r_{v},c_{v})\in G$, and $r_{v}\in[1,R]$ and $c_{v}\in[1,C]$. We refer
to the scenario as a partial-grid scenario (respectively, full-grid scenario)
when $H\subset G$ (respectively, $H=G$). Finally, given $H\subseteq G$, we
define the DP (median) $u^{*}$ as:

$u^{*}=\operatorname*{arg\,min}_{u=(r_{u},c_{u})\in G}\mathcal{C}(H,u)$ (6)

##  4 Solving MEMP with Full-Grid Scenario

In this section we give properties and then devise algorithms for solving MEMP
for the full-grid scenario. This full-grid scenario is justified by the fact
that a delivery company has to consider all the grid’s locations as potential
customers. For example, in the case of a quarantined area, as with COVID-19,
the drone can be used to deliver goods of primary necessity to all the
residents in the area. Therefore, the objective would be to find the optimal
location to set the DP in order to minimize the travel distances between any
customer’s locations and the DP.

In the following, we first discuss how to optimally solve MEMP with a full-
grid (see Section 4.2), and then we propose an approximation algorithm,
CMALL-F, that operates as if the grid is a full Manhattan grid (i.e., $K=0$)
and provides a guaranteed approximation bound of $\sqrt{2}$ (see Section
4.3.1). Then we propose two heuristics: (1) CEMB-F which assumes that all the
Manhattan destinations move to the border $B$, i.e., the grid is a fully
Euclidean grid (see Section 4.3.2), and (2) CMEB-F which does the opposite,
i.e., the grid is a fully Manhattan grid (see Section 4.3.3). In Table 1 we
compare the presented algorithms that solve MEMP in the full-grid scenario
evaluating their time complexities and guaranteed approximation bounds.

Table 1: Comparison between the algorithms that solve MEMP in the full-grid scenario. Point | Algorithm | Section | Time complexity | Approximation ratio
---|---|---|---|---
$u^{*}$ | OPT-F | 4.2 | $\mathcal{O}(\log K)$ | $1$
$u_{M}$ | CMALL-F | 4.3.1 | $\mathcal{O}(1)$ | $\sqrt{2}$
$u_{\hat{C}}$ | CEMB-F | 4.3.2 | $\mathcal{O}(1)$ | $-$
$u_{\hat{M}}$ | CMEB-F | 4.3.3 | $\mathcal{O}(1)$ | $-$

###  4.1 Properties

In the following we prove some properties that we will exploit to devise our
optimal algorithm. We first note that with a full Manhattan grid (i.e., $K=1$)
or a full Euclidean grid (i.e., $K=C$), MEMP can be trivially solved. In the
former case, MEMP has the Manhattan-median in
$u^{*}=(\overline{R},\overline{C})$ (see [44]). Note that, the median is not
unique when the values $C$ and $R$ are even. In the latter case, by using
symmetry arguments, it can be proven that MEMP has the Euclidean-median in
$u^{*}=(\overline{R},\overline{C})$. For the general case when $1<K<C$, we can
derive properties to narrow down the set of possible median point candidates.

First, we observe that the median always belongs to the middle row
$\overline{R}$ of $G$.

######  Theorem 1.

Given an EM-grid $G=(R,C,K)$, the median $u^{*}=(r^{*},c^{*})$ satisfies
$r^{*}=\overline{R}$.

###### Proof.

Let the row-cost $\Gamma(u,h)$ be the distance traversed by a drone with DP
$u=(x,y)$ to serve all the vertices on a row at distance $h$ from $u$. Note
that the function $\Gamma$ depends only on the relative distance between the
$x$-coordinate of $u$ and the row considered. There are potentially two rows
at distance $h$ from $u$: one above $u$ at row $x+h$, and one below $u$ at row
$x-h$. The crucial observation is that the two rows have exactly the same row
cost when served by $u$. For a fixed $u$, $\Gamma(u,h)$ increases with $h$. In
other words, for $h_{2}>h_{1}$, it holds that
$\Gamma(u,h_{2})-\Gamma(u,h_{1})>0$.

Now we can show that $u^{*}$ belongs to row $\overline{R}$, proving that for a
given $u=(\overline{R},j)$ and $v=(\ell,j)$, with $\ell\neq\overline{R}$, it
holds $\mathcal{C}(G,u)\leq\mathcal{C}(G,v)$. First, we consider
$\ell>\overline{R}$. In this case:

$\displaystyle\mathcal{C}(G,u)=2\sum_{x=1}^{\overline{R}-1}\Gamma(u,x)+\Gamma(u,0)$
(7)
$\displaystyle\mathcal{C}(G,v)=\sum_{x=1}^{\ell-1}\Gamma(v,x)+\sum_{x=1}^{n-\ell}\Gamma(v,x)+\Gamma(v,0)$
(8) Subtracting $\mathcal{C}(G,u)$ from $\mathcal{C}(G,v)$, one has:
$\displaystyle\sum_{x=\overline{R}}^{\ell-1}\Gamma(u,x)-\sum_{z=n-\ell+1}^{\overline{R}-1}\Gamma(u,z)\geq
0$ (9)

In Eq. (9) $x\geq\overline{R}$, while $z<\overline{R}$, and accordingly it
holds $\Gamma(u,x)\geq\Gamma(u,z)$. Hence
$\mathcal{C}(G,u)\leq\mathcal{C}(G,v)$ for any $v$. Similarly, the result can
be proven when $\ell<\overline{R}$. ∎

Recall the notion of column-cost defined in Section 3.2. Recall, also, that
$n=|G|$. Algebraically the following properties can be proven about the
column-cost:

######  Lemma 2.

  1. 1.

Both $\Delta_{E}({j})$ and $\Delta_{M}({j})$ increase with $j$;

  2. 2.

$\Delta_{M}({j})-\Delta_{M}({j-t})=t\cdot n$ for $j\geq t$;

  3. 3.

$\Delta_{E}({j})-\Delta_{E}({j-t})$ is positive and strictly monotone
increasing with $j$ for all $j>t$;

  4. 4.

$\Delta_{E}({j})-\Delta_{E}({j-t})\geq
t(\Delta_{E}({j-t+1})-\Delta_{E}({j-t}))$,
$\Delta_{E}({j})-\Delta_{E}({j-t})\leq t(\Delta_{E}({j})-\Delta_{E}({j-1}))$;

  5. 5.

$\Delta_{E}({j})<\Delta_{M}({j})\leq\sqrt{2}\Delta_{E}({j})$.

###### Proof.

We prove each point separately:

  1. 1.

Let $j_{1}>j_{2}$. From Eq. (2), it holds
$\sqrt{i^{2}+j_{1}^{2}}>\sqrt{i^{2}+j_{2}^{2}}$. From Eq. (3), it holds
$i+j_{1}>i+j_{2}$.

  2. 2.

From Eq. (3), it holds $t+2(\overline{R}-1)t=tn$.

  3. 3.

From Eq. (2), it holds $\sqrt{i^{2}+j^{2}}>\sqrt{i^{2}+(j-t)^{2}}$. Now,
differentiating with respect to $j$ one obtains:

${1}/{\sqrt{\left({i}/{j}\right)^{2}+1}}-{1}/{\sqrt{\left({i}/{j-t}\right)^{2}+1}}.$

Since $t>0$, $j>j-t$, and hence $\frac{i}{j-t}>\frac{i}{j}$, confirming the
strictly monotone increase for $j>t$.

  4. 4.

Observe that $\Delta_{E}({j})-\Delta_{E}({j-t})$ can be rewritten as the sum
of the distance of consecutive columns as
$\sum_{z=0}^{t-1}\Delta_{E}({j-z})-\Delta_{E}({j-z-1})$. By applying Property
3 with $t=1$, it holds $\Delta_{E}({j})-\Delta_{E}({j-t})\geq
t(\Delta_{E}({j-t+1})-\Delta_{E}({j-t}))$ because
$\Delta_{E}({j-z})-\Delta_{E}({j-z-1})$ is increasing with $j-z$. Similarly,
it follows: $\Delta_{E}({j})-\Delta_{E}({j-t})\leq
t(\Delta_{E}({j})-\Delta_{E}({j-1}))$.

  5. 5.

Recalling the well-known Cauchy-Schwarz inequality
$\sqrt{a^{2}+b^{2}}<(a+b)\leq\sqrt{2}\sqrt{a^{2}+b^{2}}$, it holds:
$\Delta_{E}({j})<\Delta_{M}({j})\leq\sqrt{2}\Delta_{E}({j})$. ∎

Having established in Theorem 1 that $u^{*}$ is on the row $\overline{R}$, the
potential candidates for the median are vertices $(\overline{R},c)$. For a
vertex that belongs to the middle row, say $(\overline{R},c)$, let
${\mathcal{\overline{C}}}(c)=\mathcal{C}(G,c)$. Selecting an arbitrary vertex
$u=(\overline{R},c_{u})$ as the candidate median, we can exploit the column-
cost definition and decompose the cost ${\mathcal{\overline{C}}}(c_{u})$ into
$C_{1},\ldots,C_{4}$ as defined in Eqs. (10a) and (10b). Both equations
coincide when $c_{u}=K$ because $\Delta_{M}({0})=\Delta_{E}({0})$ and
$\Delta_{M}({0})+iR=\Delta_{M}({i})$.

$\displaystyle{\mathcal{\overline{C}}}(c_{u})$
$\displaystyle=\overbrace{\underbrace{\sum\limits_{j=0}^{c_{u}-1}\Delta_{E}({j})+\sum\limits_{j=1}^{K-c_{u}}\Delta_{E}({j})}_{C_{1}=\text{cost($E$)}}+\underbrace{(C-K)\Delta_{E}({K-c_{u}})+R\sum\limits_{j=1}^{C-K}j}_{C_{2}=\text{cost($M$)}}}^{\text{if}~{}c_{u}\in
E}$ (10a) $\displaystyle{\mathcal{\overline{C}}}(c_{u})$
$\displaystyle\overbrace{\underbrace{\sum\limits_{j=1}^{K}\Delta_{E}({j})+R(K-1)(c_{u}-K)}_{C_{3}=\text{cost($E-B$)}}+\underbrace{\sum\limits_{j=0}^{c_{u}-K}\Delta_{M}({j})+\sum\limits_{j=1}^{C-c_{u}}\Delta_{M}({j})}_{C_{4}=\text{cost($M\cup
B$)}}}^{\text{if}~{}c_{u}\in M\cup B}$ (10b)

We now prove some technical results to help to further reduce the set of
median candidates. We first show, in Lemma 3, that the median cannot be “too
close” to the left border of $G$.

######  Lemma 3.

The column $c^{*}$ of the DP $u^{*}=(\overline{R},c^{*})$ of $G=(R,C,K)$
cannot be in the interval $[1,\ldots,\overline{K}-1]$, where
$\overline{K}=\lceil\frac{K}{2}\rceil$ is the column that halves the Euclidean
sub-grid.

###### Proof.

This is equivalent to saying that, if $c_{u}\in[1,\overline{K}-1]$, then
${\mathcal{\overline{C}}}(c_{u})>{\mathcal{\overline{C}}}(\left\lfloor\frac{K}{2}\right\rfloor)$.

From Eq. (10a), we notice that the cost $C_{1}$ increases if
$c_{u}<\overline{K}$ because $c_{u}$ is a sub-optimal solution for the median
of the Euclidean sub-grid in $G$. Namely, $\overline{K}$ is the median of an
EM-grid $G^{\prime}=(R,K,K)$. Moreover, by Lemma 2 (Property 3), the cost
$(C-K)\Delta_{E}({K-c_{u}})>(C-K)\Delta_{E}({\overline{K}})$ because
$K-c_{u}>\overline{K}$. ∎

Next, we show, in Lemma 4 that ${\mathcal{\overline{C}}}(c_{u})$ is convex
when $c_{u}$ varies from $\overline{K}$ to $K$.

######  Lemma 4.

Varying $c_{u}$ from $\overline{K}$ to $K$, the delivery cost function
${\mathcal{\overline{C}}}(c_{u})$ has a single minimum.

###### Proof.

Let $t$, with $\overline{K}\leq t\leq K-1$, be the vertex where the delivery
cost assumes the first minimum. Therefore from Eq. (10a) we have that:

$\displaystyle{\mathcal{\overline{C}}}(t+1)-{\mathcal{\overline{C}}}(t)$
$\displaystyle=\Delta_{E}({t+1})-\Delta_{E}({K-t})$
$\displaystyle+(C-K)(\Delta_{E}({K-(t+1)})-\Delta_{E}({K-t}))\geq 0$ (11)

To prove that the cost function has exactly one minimum in the interval
$[\overline{K},\ldots,K]$, it is sufficient to show that
${\mathcal{\overline{C}}}(t+2)-{\mathcal{\overline{C}}}(t+1)\geq 0$ given that
${\mathcal{\overline{C}}}(t+1)-{\mathcal{\overline{C}}}(t)\geq 0$.

First, observe by Lemma 2 (Property 3) that:

$\Delta_{E}({K-(t+2)})-\Delta_{E}({K-(t+1)})\geq\Delta_{E}({K-(t+1)})-\Delta_{E}({K-t})$
(12)

Then,

$\displaystyle{\mathcal{\overline{C}}}(t+2)-{\mathcal{\overline{C}}}(t+1)$
$\displaystyle=\Delta_{E}({t+2})-\Delta_{E}({K-(t+1)})$
$\displaystyle+(C-K)(\Delta_{E}({K-(t+2)})-\Delta_{E}({K-(t+1)}))$
$\displaystyle\geq\Delta_{E}({t+2})-\Delta_{E}({K-(t+1)})$
$\displaystyle+(C-K)(\Delta_{E}({K-(t+1)})-\Delta_{E}({K-t}))$
$\displaystyle=\underbrace{\Delta_{E}({t+2})-\Delta_{E}({t+1})}_{\geq
0}+\underbrace{\Delta_{E}({K-t})-\Delta_{E}({K-(t+1)})}_{\geq 0}$ (13)
$\displaystyle+{\mathcal{\overline{C}}}(t+1)-{\mathcal{\overline{C}}}(t)\geq{\mathcal{\overline{C}}}(t+1)-{\mathcal{\overline{C}}}(t)\geq
0$ (14)

∎

So, we know now that the optimal column cannot be in the interval
$[1,\ldots,\overline{K}-1]$, and that ${\mathcal{\overline{C}}}(c_{u})$ has a
single minimum in the interval $[\overline{K},K]$. We can have two cases here:
$K\leq\overline{C}$ and $K>\overline{C}$. Lemma 5 shows that when
$K\leq\overline{C}$, there is only one vertex candidate as the median on the
Manhattan side, while Lemma 6 establishes the fact that when $K>\overline{C}$,
the median cannot be in $[\overline{C},C]$.

######  Lemma 5.

Let $K\leq\overline{C}$. Varying $c_{u}$ in the Manhattan side, $K\leq
c_{u}\leq C$, the delivery cost function ${\mathcal{\overline{C}}}(c_{u})$ has
a single minimum in $\overline{C}$.

###### Proof.

Since the candidate median is on the Manhattan side, the delivery cost is
expressed by Eq. (10b). Moving the candidate from $(\overline{R},c_{u})$ to
$(\overline{R},c_{u}+1)$, the vertices on the left of $c_{u}$ (including
$c_{u}$ itself) increase by exactly one their distance from the median.
Whereas, the vertices on the right of the column $c_{u}$ decrease their
distance from the median by one. Hence, the delivery cost function decreases
as long as $K\leq c_{u}\leq\overline{C}$. ∎

######  Lemma 6.

Let $K>\overline{C}$. The value of the median candidate column cannot be
greater than $\overline{C}$.

###### Proof.

We can first exclude any candidate with $c_{u}\geq K$ because, as we seen in
Lemma 5, the delivery cost function ${\mathcal{\overline{C}}}(c_{u})$ is
increasing for $c_{u}\geq K$. Moreover, in order to exclude the candidates in
the range $[\overline{C},K]$, we first prove that
${\mathcal{\overline{C}}}(\overline{C}+1)>{\mathcal{\overline{C}}}(\overline{C})$.
Namely,

$\displaystyle{\mathcal{\overline{C}}}(\overline{C}+1)-{\mathcal{\overline{C}}}(\overline{C})$
$\displaystyle=\Delta_{E}({\overline{C}+1})-\Delta_{E}({K-\overline{C}})$
$\displaystyle-(C-1-K)\left(\Delta_{E}({K-\overline{C}})-\Delta_{E}({K-(\overline{C}+1)})\right)\geq
0$ (15)

because, by Property 4) and 3) of Lemma 2

$\displaystyle\Delta_{E}({\overline{C}+1})-\Delta_{E}({K-\overline{C}})$
$\displaystyle>\left(2\overline{C}+1-K\right)\left(\Delta_{E}({K-\overline{C}+1})-\Delta_{E}({K-\overline{C}})\right)$
$\displaystyle>(2\overline{C}+1-K)\left(\Delta_{E}({K-\overline{C}})-\Delta_{E}({K-\overline{C}-1})\right)$
$\displaystyle\geq(C-1-K)\left(\Delta_{E}({K-\overline{C}})-\Delta_{E}({K-(\overline{C}+1)})\right)$
(16)

Since we have proven in Lemma 4 that if the cost function is increasing in one
vertex of row $\overline{R}$, then it is increasing in all the vertices on its
right, there are no candidates for the median greater than $\overline{C}$. ∎

Then, Theorem 2 follows immediately from Theorem 1, and Lemmas 3, 5, and 6.

######  Theorem 2.

Suppose $u^{*}=(r^{*},c^{*})$ be the median for $G=(R,C,K)$. Let
$c,\overline{K}\leq c\leq K$ be such that

${\mathcal{\overline{C}}}(c)=\min_{i\in[\overline{K},K]}\\{{\mathcal{\overline{C}}}(i)\\}$

and, let $c^{\prime},\overline{K}\leq c^{\prime}\leq\overline{C}$ be such that

${\mathcal{\overline{C}}}(c^{\prime})=\min_{i\in[\overline{K},\overline{C}]}\\{{\mathcal{\overline{C}}}(i)\\}$

Then, $r^{*}=\overline{R}$, and

$c^{*}=\begin{cases}c&\mbox{if }K\leq\overline{C}$ \mbox{ and }
${\mathcal{\overline{C}}}(c)<{\mathcal{\overline{C}}}(\overline{C})\\\
\overline{C}&\mbox{if }K\leq\overline{C}$ \mbox{ and }
${\mathcal{\overline{C}}}(c)\geq{\mathcal{\overline{C}}}(\overline{C})\\\
c^{\prime}&otherwise\end{cases}$

∎

###  4.2 The Optimal Algorithm OPT-F

Theorem 2 directly translates to algorithm OPT-F (see the pseudo-code in
Algorithm 1) that optimally solves MEMP in the full-grid scenario.

1: if $K\leq\overline{C}$ then

12: $c\leftarrow\texttt{find-minimum}(\overline{K},K)$;

3: if ${\mathcal{\overline{C}}}(c)<{\mathcal{\overline{C}}}(\overline{C})$
then

4:

5: return $(\overline{R},c)$ {see Lemma 4}

6: else

7:

8: return $(\overline{R},\overline{C})$ {see Lemma 5}

9: end if

10: else

11: $c\leftarrow\texttt{find-minimum}(\overline{K},\overline{C})$ {see Lemma
6}

12:

13: return $(\overline{R},c)$

14: end if

Algorithm 1 The OPT-F Algorithm

Given a closed interval of column numbers, the procedure find-minimum returns
the column number, $c$, in that interval such that the total delivery cost
from the DP $(\overline{R},c)$ is the least over all the vertices
$(\overline{R},i)$ for $i$ in the given closed interval.

####  4.2.1 Time Complexity of Algorithm OPT-F

Note that the minimum in $[\overline{K},K]$, returned by invoking the
procedure find-minimum, can be found by applying a binary search due to the
unimodality proven in Lemma 4. Similarly, when $K>\overline{C}$ (Line 10), the
minimum is in the interval $[\overline{K},\overline{C}]$, i.e., a sub-interval
of $[\overline{K},K]$. Thus, as above, it can be found through the find-
minimum procedure in Line 11. The time complexity of the find-minimum
procedure is logarithmic in the width of the sub-interval where the minimum
resides. The interval has a width of $\frac{K}{2}$ in Line 1, and a width of
$\frac{C}{2}-\frac{K}{2}\leq K-\frac{K}{2}=\frac{K}{2}$ in Line 11 since
$K\geq\frac{C}{2}$. Thus, we can conclude that, in each case, the time
complexity of find-minimum is $\mathcal{O}(\log K)$.

With regards to the time complexity of the OPT-F algorithm, we observe that
for a fixed $c_{u}$, the delivery cost ${\mathcal{\overline{C}}}(c_{u})$ can
be computed in $\mathcal{O}(1)$ time by applying Eq. (10a) if the prefix sums
of the column-cost are computed in a pre-processing phase. The column-costs
$\Delta_{E}({j})$ and their prefix sums $\sum_{t=1}^{j}\Delta_{E}({t})$, for
$1\leq j\leq K$, can be computed and memorized in a vector in
$\mathcal{O}(RK+K)$ time. For $1\leq j\leq K$, the $K$ column-costs
$\Delta_{M}({j})$ and their prefix sums $\sum_{t=1}^{j}\Delta_{M}({t})$ can be
computed and memorized in a vector in $\mathcal{O}(K)$ time by using the
closed-form in Eq. (3). Then, assuming that the prefix-sums of the column-
costs are given as input to the algorithm, i.e., they are computed in a pre-
processing phase, each ${\mathcal{\overline{C}}}(c_{u})$ can be computed in
constant time. Thus, the optimal point $u^{*}$ can be computed by the OPT-F
algorithm in $\mathcal{O}(\log K)$ time.

###  4.3 Approximation and Heuristic Algorithms

In the previous section, we presented an algorithm that optimally solves MEMP
in the full-grid scenario, taking logarithmic time in the number of columns of
the grid. In this section, we present three constant-time algorithms, CMALL-F,
CEMB-F, and CMEB-F, for solving MEMP in the full-grid scenario. We also
establish an upper bound on the approximation ratio for algorithm CMALL-F. In
Section 6.2.1, we present comparative empirical performance evaluation of
these three approximation algorithms using synthetic data.

####  4.3.1 Algorithm CMALL-F

Algorithm CMALL-F returns the Manhattan-median
$u_{M}=(\overline{R},\overline{C})$ as the DP. That is, CMALL-F operates as if
the grid is a full Manhattan grid (i.e., $K=1$) This algorithm is sub-optimal
for $1<K<C$ providing a guaranteed approximation bound of $\sqrt{2}$, while it
is optimal when $K=1$ or $K=C$.

######  Lemma 7.

The CMALL-F algorithm provides a $\sqrt{2}$ approximation ratio when $1<K<C$.

###### Proof.

Let $u^{*}$ be the median for $G=(R,C,K)$. Since the Euclidean column-cost
$\Delta_{E}({j})$ is smaller than the Manhattan column-cost $\Delta_{M}({j})$,
$\mathcal{C}(G,u_{M})\leq\mathcal{C}_{M}(G,u_{M})$, where
$\mathcal{C}_{M}(G,u_{M})$ is the cost when $K=1$. Moreover,
$\mathcal{C}(G,u^{*})>\mathcal{C}_{E}(G,u^{*})$ because the Manhattan distance
is at least as much as the Euclidean distance. Then,
$\mathcal{C}_{E}(G,u^{*})>\mathcal{C}_{E}(G,u_{M})$ because $u_{M}$ is the
Euclidean-median of $G=(R,C,C)$. Thus, by the Cauchy-Schwarz inequality in
Lemma 2:

$\frac{\mathcal{C}(G,u_{M})}{\mathcal{C}(G,u^{*})}<\frac{\mathcal{C}_{M}(G,u_{M})}{\mathcal{C}_{E}(G,u_{M})}\leq\sqrt{2}.$

The CMALL-F algorithm finds the point $u_{M}$ in constant time, so its time
complexity is $\mathcal{O}(1)$.

####  4.3.2 Algorithm CEMB-F

Algorithm CEMB-F solves MEMP with full-grid, and selects the DP
$u_{\hat{C}}=(\overline{R},\mu)$ with

$\mu=\frac{(\sum_{i=1}^{K}i)+(C-K)K}{C}=\frac{K(2C-K+1)}{2C}$ (17)

CEMB-F imagines that all the Manhattan destinations move on the border $B$.
Thus, the grid becomes a full Euclidean grid $E^{\prime}$ with $K$ columns,
whose rightmost column has multiplicity $w_{v}=C-K$. So, minimizing Eq. (10a)
is the same as finding the Euclidean-median of $E^{\prime}$. Unfortunately,
there is no closed form to compute the exact Euclidean-median for a set of
points111We only know the Euclidean-median of a Euclidean grid $G=(R,C,C)$,
whose columns have all multiplicity $1$..

Algorithm CEMB-F finds the point $u_{\hat{C}}$ using a constant number of
operations.

####  4.3.3 Algorithm CMEB-F

Algorithm CMEB-F imagines that all the Euclidean destinations are moved on the
border $B$. Thus, the grid becomes a full Manhattan grid $M^{\prime}$ with
$C-K$ columns, whose leftmost column has multiplicity $w_{v}=K$. Therefore,
minimizing Eq. (10b) is the same as finding the Manhattan-median of the grid
$M^{\prime}$. Algorithm CMEB-F solves MEMP with full-grid, and selects the DP
$u_{\hat{M}}=(\overline{R},\mu)$ with $\mu=K$ if $K\geq\overline{C}$ or with
$\mu=\overline{C}$ if $K<\overline{C}$. In either case, $u_{\hat{M}}$ is
computed in constant time.

##  5 Solving MEMP with Partial-Grid Scenario

In this section, we focus on MEMP in the partial-grid scenario. In this case,
we are given $H\subset G$ as the $n$ delivery points, i.e., $|H|=n$.

We first discuss the trivial cases where $K=1$, i.e., a full Manhattan grid,
and $K=C$, i.e., a full Euclidean grid (see Section 5.1). The rest of the
sections deals with the general cases.

In Sections 5.2 and 5.3 we propose two heuristic algorithms, CEMB-P and
CMEB-P, assuming the optimal DP is in the Euclidean side of the grid, and in
the Manhattan side of the grid, respectively. Note that neither of these two
algorithms returns the optimal DP. Algorithm CEMB-P returns the best DP
candidate in $E$, whereas CMEB-P returns the best DP candidate in $M$.

Then, in Section 5.4 we combine the above two heuristic algorithms and present
algorithm OPT-P to find the optimal median, the DP, in $G$.

In Sections 5.5 and 5.6 we present two algorithms, S-OPT-P and CMALL-P,
respectively, that each compute a sub-optimal DP, but are more efficient than
Algorithm OPT-P.

To summarize the above, we tabulate in Table 2 the presented algorithms that
solve MEMP in the partial-grid scenario along with their time complexities and
guaranteed approximation bounds. In Section 6.2.2, we present comparative
empirical performance evaluation of these algorithms using synthetic data,
while in Section 6.3, we compare these algorithms using quasi-real data.

Table 2: Comparison between the algorithms that solve MEMP in the partial-grid scenario. Point | Algorithm | Section | Time complexity | Approximation ratio
---|---|---|---|---
$u_{\hat{C}}$ | CEMB-P | 5.2 | $\mathcal{O}(nR\log K)$ | $-$
$u_{\hat{M}}$ | CMEB-P | 5.3 | $\mathcal{O}(nR)$ | $-$
$u^{*}$ | OPT-P | 5.4 | $\mathcal{O}(nR\log K)$ | $1$
$\overline{u}$ | S-OPT-P | 5.5 | $\mathcal{O}(n\log R\log K)$ | $-$
$u_{M}$ | CMALL-P | 5.6 | $\mathcal{O}(n)$ | $-$

###  5.1 $K=1$ or $K=C$

We first discuss how to optimally solve MEMP to serve a subset of customers on
a grid that is not mixed, i.e., either Manhattan or Euclidean. In a Manhattan
grid (i.e., an EM-grid with $K=1$), given a subset $H\subset G$ of customers,
MEMP with partial-grid scenario can be trivially solved in $\mathcal{O}(|H|)$
time by selecting $u^{*}=(r_{u^{*}},c_{u^{*}})$ where $r_{u^{*}}$ is the
median of the row coordinate of the customers and $c_{u^{*}}$ is the median of
the column coordinate of the customers (see [44]).

In the literature, there are many results concerning the problem, called
geometric median, of determining the DP that minimizes the sum of distances
between the points of a given set $H\subset\mathbb{R}^{d}$ and the DP itself.
The geometric median is denoted as the Euclidean-median in two-dimensional
space. Although the Euclidean-median is unique and the sum $\mathcal{C}(H,u)$
of the distances from each customer to the DP $u$ is positive and strictly
convex in $\mathbb{R}^{d}$, as proven by [45], there is no exact and closed
expression for the Euclidean-median of an arbitrary set $H$ of real points. In
a Euclidean grid (i.e., an EM-grid with $K=C$), MEMP with partial-grid
scenario can be solved more easily because the candidate positions in the
plane are just the vertices of $G$, and the number of attempts to determine
the single Euclidean-median is limited by the fact that EM-grid is formed by
$R$ rows and $C$ columns. So, a trivial solution for MEMP with partial-grid
takes $\mathcal{O}(RC|H|)$ time because for each position $u$ of the grid the
cost $\mathcal{C}(H,u)$ can be computed in time $\mathcal{O}(|H|)$. Exploiting
the fact that, in $\mathbb{R}^{2}$, $\mathcal{C}(H,u)$ is positive and
strictly convex when $u$ moves on a single row of the grid, the position in a
row that provides the minimum cost $\mathcal{C}(H,u)$ can be computed by a
slightly modified binary search. So, MEMP with a subset of customers on a
Euclidean grid can be solved in $\mathcal{O}(R\log C|H|)$ time because for
each row of the grid only $\log C$ candidates are tested.

In the next section, for the general case, i.e., the mixed-grid with $1<K<C$,
we leverage the observation that the median $u^{*}$ resides either in $M$ or
in $E$.

###  5.2 The CEMB-P Algorithm

As noted above, CEMB-P solves MEMP assuming that $u^{*}\in E$ and returns
$u_{\hat{C}}\in E$. Algorithm CEMB-P selects as the DP the point in the
Euclidean side $E$ that returns the minimum cost.

For any DP $u\in E$, by Lemma 1, for each point $v$ in $H_{M}$ (the subset of
$H$ that contains the vertices in the Manhattan side), the drone must travel
horizontally from the border $B$ to $v$. So, the drone covers the same fixed
cost in the Manhattan grid regardless of the position of $u\in E$. This
motivates our algorithm whose pseudo-code is presented in Algorithm 2.

Similar to CEMB-F for the full-grid scenario, we construct the multi-set
$H_{M}^{\prime}$ consisting of the projections of all the points in $H_{M}$ to
the border $B$ (Algorithm 2, Line 2). Then, we compute the point
$u_{\hat{C}}\in E$ that minimizes the sum of the distances to the points in
the multi-set $H^{\prime}=H_{E}\cup H_{M}^{\prime}$, taking into account the
multiplicity of points in $H^{\prime}$. Towards this, for each row (Algorithm
2, Line 3) we evaluate the minimum from the first to the $K^{th}$ column, and
update the overall minimum, if necessary. By exploiting the fact that in the
Euclidean space the function $\mathcal{C}(H,u)$ to minimize is positive and
strictly convex (see [45]), we know that $\mathcal{C}(H,u)$ is unimodal when
fixing a row and varying the columns. Hence, we can calculate the minimum on
each row by performing a time-efficient binary search that requires
logarithmic time (Algorithm 2, Line 2). Since the cost of movement in $M$ does
not change the minimum, the returned value, $u_{\hat{C}}$, is the best DP in
$E$.

11: $H_{M}^{\prime}\leftarrow\\{(r_{u},K)\in B\mid u\in
H_{M}\\},H^{\prime}\leftarrow H_{E}\cup H_{M}^{\prime}$;

22: $u_{\hat{C}}\leftarrow\varnothing,cost\leftarrow+\infty$;

3: for $i\in 1,\ldots,R$ do

34: $u^{*}_{i}\leftarrow\texttt{find-minimum-on-row}(i,H^{\prime},1,K)$;

5: if $\mathcal{C}(H^{\prime},u^{*}_{i})<cost$ then

46: $u_{\hat{C}}\leftarrow
u^{*}_{i},cost\leftarrow\mathcal{C}(H^{\prime},u^{*}_{i})$;

7: end if

8: end for

9:

10: return $u_{\hat{C}}$

Algorithm 2 The CEMB-P Algorithm

About the time complexity, due to the fact that there are $n$ points to serve,
and since we perform $R$ binary searches (one for each row of a Euclidean grid
with $K$ columns), the total cost of CEMB-P is $\mathcal{O}(nR\log K)$.

###  5.3 The CMEB-P Algorithm

Algorithm CMEB-P solves MEMP assuming that $u^{*}\in M$ and returns
$u_{\hat{M}}\in M$.

For any DP, $u$ in $M$, by Lemma 1, the drone has to fly through the
projection of $u\in B$ to serve any point of $H_{E}$. Thus, for any DP $u$,
there is one intermediate point in $B$ for all drone paths from $u$ to points
in $H_{E}$. Note that the intermediate point does not depend on the column
$c_{u}$ of the DP $u$. This motivates our algorithm whose pseudo-code is
presented in Algorithm 3.

We first compute the column median, $\chi$, of the points in $H^{\prime}$,
consisting of the points in $H_{M}$ and the points in $H_{E}$ moved to the
border. The function column-median takes into account the possible
multiplicity of points in $H^{\prime}$. Since we are only concerned with the
column numbers, we can move each point in $H_{E}$ to any point on the border;
we move all of them to $(1,K)$ (Algorithm 3, Lines 3 and 3).

Since the median in $M$ can reside in $R$ different rows, we have $R$
candidate intermediate points, which are the points of the border $B$. For
each row $i$, we construct the multi-set $H_{E}^{\prime}(i)$ consisting of all
the points in $H_{E}$ moved to the intermediate point $(i,K)$ on the border
$B$ (Algorithm 3, Line 3). We then consider the point $u_{i}^{*}=(i,\chi)\in
M$ as the median of the points in the multi-set $H^{\prime}=H_{M}\cup
H_{E}^{\prime}(i)$. With $u_{i}^{*}$ as the DP, the cost of delivery is the
sum of two costs: (1) the cost of flying the drone between $u_{i}^{*}$ and
each point in $H^{\prime}$, with all distances calculated according to the
Manhattan metric, and (2) the cost of flying the drone between each point in
$H_{E}$ and the intermediate point corresponding to $u_{i}^{*}$, with all
distances calculated according to the Euclidean metric (Algorithm 3, Line 8).

The algorithm returns as $u_{\hat{M}}$ the intermediate point that witnesses
the least cost, over all the possible intermediate points, and is thus the
best DP in $M$.

11: $u_{\hat{M}}\leftarrow\varnothing,c\leftarrow+\infty$;

22: $H_{E}^{\prime}(1)\leftarrow\\{(r_{u},1)\in B\mid u\in H_{E}\\}$,
$H^{\prime}\leftarrow H_{E}^{\prime}(1)\cup H_{M}$;

33: $\chi\leftarrow\texttt{column-median}(H^{\prime})$;

4: for $i\in 1,\ldots,R$ do

45: $H_{E}^{\prime}(i)\leftarrow\\{(r_{u},i)\in B\mid u\in H_{E}\\}$,
$H^{\prime}\leftarrow H_{E}^{\prime}(i)\cup H_{M}$;

56: $u^{*}_{i}\leftarrow(i,\chi)$;

67:
$C_{i}\leftarrow\mathcal{C}(H^{\prime},u^{*}_{i})+cost(H_{E}\rightarrow(i,K))$;

8: if $C_{i}<c$ then

79: $u_{\hat{M}}\leftarrow u^{*}_{i}$;

810: $c\leftarrow C_{i}$;

11: end if

12: end for

13:

14: return $u_{\hat{M}}$

Algorithm 3 The CMEB-P Algorithm

As for the complexity of the algorithm, we compute the column median of $n$
points once. Then, for every row $i$, the cost $cost(H_{E}\rightarrow(i,K))$
has to be computed, and this requires $\mathcal{O}(|H_{E}|)$ time. Thus, the
algorithm’s complexity is $\mathcal{O}(n+R|H_{E}|)$.

###  5.4 The Optimal OPT-P Algorithm

Having computed the best DPs on both sides of EM-grid we can now optimally
solve MEMP with partial-grid scenario.

The OPT-P algorithm finds the optimal point $u^{*}$ by comparing the best
point between $u_{\hat{C}}$ and $u_{\hat{M}}$. Given that Algorithm CEMB-P
returns the best DP in $E$, and Algorithm CMEB-P returns the best DP in $M$,
the simple idea of OPT-P is to compare these two points and return the best
one, as follows:

$u^{*}=\operatorname*{arg\,min}\\{\mathcal{C}(H,u_{\hat{C}}),\mathcal{C}(H,u_{\hat{M}})\\}.$
(18)

About the time complexity, recalling that CEMB-P takes $\mathcal{O}(nR\log K)$
and CMEB-P takes $\mathcal{O}(nR)$, the overall time complexity of OPT-P is
$\mathcal{O}(nR\log K)$.

###  5.5 The S-OPT-P Algorithm

Algorithm S-OPT-P applies the same strategy as OPT-P, i.e., of comparing the
best among two points: the best one in the Euclidean grid, and the best one in
the Manhattan grid. Nonetheless, it uses slightly different versions of both
CEMB-P and CMEB-P.

In the modified version of Algorithm CEMB-P, we perform two binary searches –
one on the rows plus the one on the columns up to the $K^{th}$ column. In the
modified version of Algorithm CMEB-P, we perform a single binary search on the
rows. Thus, the time complexities of the two algorithms are
$\mathcal{O}(nR\log K)$ and $\mathcal{O}(nR)$, respectively.

Thus, the overall time complexity of S-OPT-P is $\mathcal{O}(n\log R\log K)$.
However, this strategy does not guarantee that the returned point
$\overline{u}$ is optimal.

###  5.6 The CMALL-P Algorithm

Essentially, Algorithm CMALL-P ignores the border that separates the Euclidean
and the Manhattan sides, and computes the DP $u_{M}$ as if $G=(R,C,1)$. The
algorithm returns

$u_{M}=(\tilde{r}_{H},\tilde{c}_{H}),$ (19)

where $\tilde{r}_{H}$ and $\tilde{c}_{H}$ are the individual medians of the
rows and the columns, respectively, of the points in $H$ (see [44]). Note
that, also in this scenario, the median is not unique if $|H|$ is even.

Although CMALL-P optimally solves MEMP in the case of $G=(R,C,1)$, it is sub-
optimal in the general case. Nonetheless, it works in linear time (see [46])
with respect to the number of points.

##  6 Performance Evaluation

In this section, we empirically compare the performance of our algorithms in
terms of the quality of the solution (i.e., delivery cost), and their running
times, for solving MEMP in both scenarios.

###  6.1 Settings and Parameters

We implemented our algorithms in Python language version 3.9, and run all the
instances on an Intel i7-860 computer with $12\mathrm{\,GB}$ of RAM. In order
to evaluate our proposed algorithms for solving MEMP in both scenarios, we
rely on synthetic and quasi-real delivery areas.

For the synthetic case (Section 6.2), we set different layouts by varying
$R,C\in\\{50,\ldots,400\\}$ and $1\leq K\leq C$. Then, we compare the
algorithms with respect to the optimal one, and we plot the experimental ratio
$\rho=\frac{\mathcal{C}(H,\tilde{u})}{\mathcal{C}(H,u^{*})}\geq 1$. In other
words, for $H\subseteq G$, $\rho$ is the ratio between the total cost for
serving all the required customers from the DP $\tilde{u}$ as returned by the
tested algorithm, and the total cost of the optimal solution where the DP used
is $u^{*}$. When testing the full-grid scenario, we compare CMALL-F, CEMB-F,
and CMEB-F with respect to the optimal algorithm OPT-F, while when testing the
partial-grid scenario, we compare CMALL-P, CEMB-P, CMEB-P, and S-OPT-P with
respect to the optimal algorithm OPT-P.

Moreover, in the partial-grid scenario, we uniformly generate $n=|H|$ random
positions inside the grid with $n=\\{5,\ldots,100\\}$, and then return the
average ratio (along with the standard deviation) on $33$ random instances.
Also in the partial-grid scenario, given a setting with $n$ random customers,
we evaluate the algorithms when balancing the quantities $n_{E}$ and $n_{M}$
with respect to a certain fraction
$p=\\{0,\frac{1}{3},\frac{1}{2},\frac{2}{3},1\\}$ on $n$, such that
$n_{E}=p\cdot n$ and $n_{M}=(1-p)\cdot n$, with $n=n_{E}+n_{M}$.

For the quasi-real case (Section 6.3), we only test the more general partial-
grid scenario in random instances taken from real cities, like the ones shown
in Figure 1. For these examples, we approximately extract the actual EM-grid
from the map, and then we run our proposed algorithms. Obviously, real cities
cannot be exactly modeled as EM-grids due to the fact that roads and buildings
can be arbitrarily made by people. However, we have found interesting examples
and attempted to perform our algorithms on these layouts. Clearly, in the
aforementioned grids of customers, some houses or skyscrapers can be missing.

###  6.2 Results with Synthetic Data

####  6.2.1 Full-Grid Scenario

We first analyze our empirical results of the performance of the algorithms
with respect to the delivery costs, and then with respect to the running
times.

Delivery Costs

Figure 3 compares the algorithms when solving MEMP with the full-grid scenario
reporting, for each plot, the ratio $\rho$ between the total cost of the
tested algorithm and the optimal total cost.

Figure 3: All the algorithms with the full-grid scenario.

CMALL-F performs very well since its DP $u_{M}$ is always very close to
$u^{*}$. On the other hand, as expected, CEMB-F and CMEB-F are heavily
affected by the value of $K$. In fact, for small values of $K$, CEMB-F reports
bad results, while when $K$ increases the ratio $\rho$ tends to $1$.

We note that the performance of CEMB-F is almost a reflection, in the vertical
line $\frac{K}{C}=0.5$, of the performance of CMEB-F. In particular, in Figure
3, we can observe that CMALL-F and CMEB-F perform similarly when
$0\leq\frac{K}{C}\leq 0.5$ (the two lines, i.e., the orange and the blue ones,
almost coincide), and CMALL-F and CEMB-F perform similarly when
$\frac{K}{C}\leq 0.5\leq 1$ (orange and green lines).

We also note that, when $R=C$ (Figure 3, first row), the worst cases of CEMB-F
and CMEB-F, i.e., the worst $\rho$ values exhibited, almost coincide; when
$R>C$ (Figure 3, second row), the worst case of CMEB-F is slightly better than
that of CEMB-F; and when $R<C$ (Figure 3, third row) the worst case of CMEB-F
is worse than that of CEMB-F.

It is interesting to note that the performance of each algorithm is better
when $R>C$, than when $R<C$.

Figure 4: CMALL-F vs BEST-F with the full-grid scenario.

The almost symmetric performance, across $K/C=0.5$, exhibited by CEMB-F and
CMEB-F motivates the design of a hybrid algorithm, BEST-F, that simply returns
the best result among the two. In Figure 4 we compare CMALL-F against the best
among CEMB-F and CMEB-F (highlighted as BEST-F). Since the experimental
results are very close to the optimum ratio, which is $1$, here we reduce the
scale along the $y$-axis. In general, BEST-F almost always outperforms
CMALL-F, and it is very close to OPT-F. For instance, when $R=150$ and $C=50$,
BEST-F is constantly $1$, and so we can state that when $R\gg C$, BEST-F is
comparable to OPT-F. With the scaled $y$-axis in Figure 4, it becomes clear
CMALL-F has the worst performance when $\frac{K}{C}\approx 0.5$. Nevertheless,
under these circumstances, the ratio $\rho$ is very low. It is worth noting
that in the case of $C\gg R$, BEST-F is not better than CMALL-F (e.g.,
$R=50,C=400$). In fact, when $C=2R$ BEST-F is preferable, but when $C=8R$ (or
more) it seems that CMALL-F is more consistent. However, in hypothetical real
cases with a reasonable mix of rows and columns, BEST-F basically performs
like the optimal OPT-F, albeit more efficiently, since BEST-F is a constant
time algorithm. Finally, it can be seen that CMALL-F is always below the
guaranteed threshold $\sqrt{2}$ on each plot.

Figure 5: Running time of all the algorithms with the full-grid scenario.

Running Times

In Figure 5 we report the experimental running time (in milliseconds) of all
the algorithms with the full-grid scenario. In particular, the constant-time
algorithms like CMALL-F take, on average, $2$–$6\mathrm{\,ms}$, and therefore
the graphs of their running times almost coincide. On the other hand, the
optimal logarithmic-time algorithm OPT-F takes much more time – in the order
of tenths of a second (in the worst case, $200\mathrm{\,ms}$). The
experimental running time of the OPT-F algorithm is low when the grid is
either almost Manhattan or almost Euclidean, i.e., when $\frac{K}{C}\to 0$ or
$\frac{K}{C}\to 1$, respectively.

In the following we will show that the behavior of OPT-F is consistent with
the claimed time complexity $\mathcal{O}(\log K)$. To better understand this
behavior we discuss the following example (for almost Manhattan grids): When
$K$ is small with respect to $C$, say $\frac{K}{C}=0.25$ (see Figure 5),
Algorithm OPT-F evaluates the column number, in the interval
$[\overline{K},K]$, with the least cost (Line 1 of Algorithm 1). Specifically,
it performs an efficient binary search in the interval $[\overline{K},K]$. If
we assume that $\frac{K}{C}=0.25$ or $\frac{K}{C}=0.5$ (Line 1 of Algorithm
1), then the interval $[\overline{K},K]$ has width $\frac{K}{2}$. Since when
$C$ is fixed, varying $\frac{K}{C}$ from $0.25$ to $0.5$, the interval width
$\frac{K}{2}$ increases from $\frac{C}{8}$ to $\frac{C}{4}$, the number of
steps required by the OPT-F algorithm increase accordingly to the logarithm of
$\overline{K}$.

When $K$ is large with respect to $C$, i.e., the grid is almost all Euclidean,
say $\frac{K}{C}=0.75$ (see Line 11 of Algorithm 1), the OPT-F algorithm
performs an efficient binary search in the interval
$[\overline{K},\overline{C}]$. Recall that $\overline{C}=\frac{C}{2}$. If we
assume that $\frac{K}{C}=0.75$, then $K=\frac{3}{4}C$ and
$\overline{K}=\frac{K}{2}=\frac{3}{8}C$, and therefore the interval
$[\overline{K},\overline{C}]$ has width $\frac{C}{8}$. The time complexity is
therefore logarithmic in $\frac{C}{8}$ as it was in the case
$\frac{K}{C}=0.25$. Thus, the time complexity is comparable with that for
$\frac{K}{C}=0.25$, as reported in Figure 5.

It is very important to note that, the plots in Figure 5 also include the pre-
processing time required by the Algorithm OPT-F; this pre-processing has time
complexity $\mathcal{O}(RK+K)$ and hence increases with $R$ and $K\leq C$. In
Figure 5, the pre-processing phase impacts more when $C=100$ because its time
complexity depends on $K$, which can be $C$ in the worst case scenario.

####  6.2.2 Partial-Grid Scenario

As before, we will first analyze our empirical results of the performance of
the algorithms with respect to the delivery costs, and then with respect to
the running times.

Delivery Costs

Figure 6 compares the algorithms when solving MEMP with the partial-grid
scenario. In particular, given the bad behavior of Algorithm CEMB-F for large
values of $K$ in the full-grid scenario (see Figure 3), for its partial-grid
version, i.e., CEMB-P, we only report the results using small values of $K$.
For a similar reason, we report the results for Algorithm CMEB-P only for
large values of $K$. Moreover, recall that the best DP returned by CEMB-P and
CMEB-P is the optimum.

Figure 6: All the algorithms when varying $n$ with the partial-grid scenario.

As before, in Figure 6 we present our empirical results in three groups, –
$R=C$, $R>C$, and $R<C$. As expected, both the CEMB-P and CMEB-P algorithms
have different behavior, and one is more suitable than the other depending on
the particular value of $K$. When $n$ is small, all the algorithms are
variable and the standard deviation is large. Much more stable results can be
observed when $n$ increases. This is due to the fact that for larger values of
$n$ the subset $H\subset G$ is “denser”, and all the algorithms are less
affected by the randomness. In general, the algorithms seem to perform better
when $R>C$. Of particular note is the behavior exhibited by Algorithm S-OPT-P–
it performs almost as well as Algorithm OPT-P! We discuss this further in
Section 6.2.3 below.

Figure 7: All the algorithms when fixing $R=50$ and $C=50$, and varying $p$
with the partial-grid scenario.

In Figure 7 we present our empirical results for all the algorithms with the
partial-grid scenario when varying the parameter $0\leq p\leq 1$ such that
$n_{E}=n\cdot p$ and $n_{M}=n\cdot(p-1)$. For simplicity, in Figure 7 we only
depict the squared layout, in which on each row we change the border position.
Not surprisingly, when $n$ is low the variability (standard deviation) is
high, and algorithms like CMALL-P return a good ratio $\rho$ (average). The
other two algorithms, i.e., CEMB-P and CMEB-P, have the usual symmetric
behavior. However, in this setting, such behavior is much more evident since
we vary the number of deliveries to distribute in the two sub-areas. So, for
small values of $p$, CMEB-P is almost optimum, and for large values of $p$
CEMB-P is almost optimum. Instead, when $p\approx 0.5$, the two algorithms
lose their efficacy. This can be seen for the three values of $K$ shown in
Figure 7. Once again, surprising, but expected based on the results presented
in Figure 6, S-OPT-P performs almost as well as the optimum algorithm.

Figure 8: Running time of all the algorithms with the partial-grid scenario.

Running Times

Finally, in Figure 8 we report the experimental running time (in milliseconds)
of all the algorithms with the partial-grid scenario. Notice that the plots
have a different scale in the $y$-axis, i.e., the ones in the first row have
$30\mathrm{\,ms}$, while the ones in the second row have $60\mathrm{\,ms}$.
This has been done to emphasize the linear dependency on $n$: when $n$
doubles, the time performance on the $y$-axis doubles too, but since the scale
of $y$ is halved, the behavior appears to be the same. As expected, the
fastest algorithm is CMALL-P whose time complexity is $\mathcal{O}(n)$,
followed by CMEB-P with a time complexity of $\mathcal{O}(nR)$ (see Table 2).
In Figure 8, it is difficult to appreciate the actual difference between
CMALL-P and CMEB-P, but, for instance, when $\frac{K}{C}=0.5$, on average,
CMALL-P takes $0.061\mathrm{\,ms}$, while CMEB-P takes $0.786\mathrm{\,ms}$
(i.e., $10\times$ more than CMALL-P).

The running time of the CMEB-P is “linear” since its cost does not depend on
the value of the border $K$. The sub-optimal algorithm S-OPT-P experimentally
performs very well (theoretically it takes $\mathcal{O}(n\log R\log K)$ in
time), in the order of $5$–$15\mathrm{\,ms}$, with respect to the number of
deliveries $n$. The CEMB-P algorithm, which has a time complexity of
$\mathcal{O}(nR\log K)$, depends on the binary search function invoked $R$
times. The binary search itself depends on the size of the border $K$, and
therefore its running time increases (by a factor of $\log K$) when the EM-
grid tends to become more “Euclidean”.

In the full-grid scenario (Figure 5), the algorithm that exploits the binary
search (i.e., OPT-F) has a parabolic trend, justified by the fact that when
$K$ is either small or large with respect to $C$, Algorithm 1 executes
differently. In this case (partial-grid scenario), though, we do not have a
similar behavior, and the time complexity of Algorithm CEMB-P is simply
affected by the factor $\log K$. Nevertheless, it is still possible to observe
a light drop for CEMB-P when $\frac{K}{C}\approx 0.75$. Obviously, the optimal
algorithm OPT-P, whose time complexity is $\mathcal{O}(nR\log K)$, takes the
most time (as seen in Figure 8) since it runs both CEMB-P and CMEB-P, and
returns the best point among the ones outputted by both.

####  6.2.3 Noteworthy Observations

Based on the theoretical and empirical analyses above, we make the following
noteworthy observations about our algorithms for solving MEMP in the partial-
grid scenario.

In general, Algorithm CMALL-P has the best trade-off with respect to obtained
performance and time complexity. In fact, when varying the number $n$ of
deliveries, on average, its ratio $\rho$ is always lower than $1.03$.
Considering its time complexity – linear with respect to the number of
deliveries – CMALL-P is a very good compromise. As noted above, both the
CEMB-P and CEMB-P algorithms have different behavior, and one is more suitable
than the other depending on the particular value of $K$. However, as proven in
Section 4.1, their combination, by returning the best among the two, ensures
finding the optimal DP. So, despite their swinging behavior and relatively
high complexity, they can be used to find the optimal DP.

We noted above the surprising, almost-optimal behavior of the sub-optimal
algorithm S-OPT-P: On the $33$ random instances, for each combination of rows,
columns, border, and number of deliveries, the returned DP was also the
optimal one. For this reason, on each plot of Figure 6, S-OPT-P exhibits a
completely flat line. However, on an extensive “brute-force” campaign with
$10,000+$ random instances, we found $17$ instances where the returned point
$\overline{u}$ was different from the optimal $u^{*}$. Nevertheless, in these
$17$ instances the worst ratio $\rho$ was $1.005$. Of course, we have found
counterexamples that show that S-OPT-P is not optimum, but it performs
extremely well with a time complexity much lower than that of OPT-P.

###  6.3 Results with Quasi-real Data

In this section, we test our algorithms in the quasi-real case. Specifically,
we evaluate the performance of our partial-grid scenario algorithms applied on
a real-world map. Clearly, in this case, we still generate random instances of
customers.

####  6.3.1 Settings

In this section, we first discuss how the delivery maps are extracted from
real cities, and then which commercial drones are capable of delivering
products in these areas.

Map Extraction

In order to evaluate our algorithms in a real-world case, we initially needed
to extract the EM-grids from real cities. Figure 9 illustrates how to roughly
perform this extraction task by modeling Chicago, New York, and Miami, in the
top, middle, and bottom, respectively.

(a) Chicago.

(b) The corresponding EM-grid $G_{1}=(8,14,6)$.

(c) New York.

(d) The corresponding EM-grid $G_{2}=(7,24,1)$.

(e) Miami.

(f) The corresponding EM-grid $G_{3}=(9,20,20)$.

Figure 9: The original cities and the corresponding extracted EM-grids.

In Figure 9a, there is the top view of a portion of the city of Chicago, while
in Figure 9b the corresponding EM-grid, i.e., $G_{1}=(8,14,6)$. The average
length of a “block” is $\approx 120\mathrm{\,m}$, and so the distance among
two adjacent vertices can be set. So, this map $G_{1}$ has a size of
$840\mathrm{\,m}\times 1560\mathrm{\,m}$. Roughly, the river from the top to
the bottom of Figure 9a splits the portion of Chicago in two contiguous areas,
i.e., the left one with relatively low buildings (although some tall buildings
are present), and the right one with lots of skyscrapers (although some low
buildings are present). The dashed line in Figure 9b delimits the border
between the two areas.

In Figure 9c, there is the top view of a portion of the city of New York,
while in Figure 9d the corresponding EM-grid, i.e., $G_{2}=(7,24,1)$. This
means that the extracted EM-grid is a full Manhattan grid. In this case,
differently from Chicago, we can observe that the length of a block on
“streets” (horizontal) is $\approx 80\mathrm{\,m}$, while the length of a
block on “avenues” (vertical) is $\approx 150\mathrm{\,m}$. So, this map
$G_{2}$ has a size of $660\mathrm{\,m}\times 2530\mathrm{\,m}$. Since in our
model both the lengths must be equal, we average the two and set a value of
$110\mathrm{\,m}$.

Finally, in Figure 9e, there is the top view of a portion of the city of
Miami, while in Figure 9f the corresponding EM-grid, i.e., $G_{3}=(9,20,20)$.
This means that the extracted EM-grid is a full Euclidean grid. Even the pure
Euclidean maps are hard to find since roads and houses can have different
sizes, lengths, and so on. The portion of Miami city reported in Figure 9e
contains some irregularities, especially in the right portion. Here, the
average distance among two adjacent houses is in the order of $\approx
25\mathrm{\,m}$. So, this map $G_{3}$ has a size of $200\mathrm{\,m}\times
475\mathrm{\,m}$.

Drone Selection

Previously, we have seen how to extract an EM-grid from a real city.
Therefore, the next step is the drone selection. There are a plethora of off-
the-shelf drones that are available at the moment, but only a very few have
the required characteristics for performing deliveries. One of the most common
drones able to do deliveries with some autonomy is the DJI Matrice 300 RTK
(briefly Matrice). According to [47], this drone can carry up to
$2.7\mathrm{\,kg}$, can fly up to $8\mathrm{\,km}$ away from the remote
controller, can fly for $30\mathrm{\,min}$ (with maximum payload), and can fly
up to $17\mathrm{\,m/s}\approx 60\mathrm{\,km/h}$ in P-mode (i.e., the
Positioning mode recommended in autonomous flights). Clearly, the temporal
bound of $30\mathrm{\,min}$ also depends on the drone’s speed and the payload
as well as the current weather conditions. Nevertheless, we still continue to
use these numbers as a reference. So, the Matrice can fly a distance of
approximately $30\mathrm{\,km}$. Recalling that the previous three grids
$G_{1}$, $G_{2}$, and $G_{3}$ have an area of approximately
$1.2\mathrm{\,km^{2}}$, $1.5\mathrm{\,km^{2}}$, and $0.1\mathrm{\,km^{2}}$,
respectively, we believe that the Matrice can operate in these neighborhoods
of the cities.

####  6.3.2 Results

Finally, Figure 10 compares the algorithms in the partial-grid scenario in a
quasi-real case. In particular, in the $x$-axis we report the number of
deliveries $n$, while in the $y$-axis we report the traveled drone’s distance
in $\mathrm{\,km}$. As aforementioned, the selected Matrice drone can fly for
approximately $30\mathrm{\,km}$, and for this reason we put a solid red line
on that threshold in each plot of Figure 10. Moreover, recall that $G_{1}$,
$G_{2}$, and $G_{3}$ correspond to the map of Chicago (mixed EM-grid), New
York (full Manhattan EM-grid), and Miami (full Euclidean EM-grid),
respectively.

Figure 10: Traveled drone’s distance when comparing all the algorithms in the
quasi-real case. The EM-grids are the ones extracted from the maps in Figure
9.

Our first observation is that the differences among the algorithms, in terms
of distance ($\mathrm{\,km}$) traveled, is small, with a couple of exceptions.
In fact, as we discussed in the synthetic case, it is not profitable at all to
perform the CEMB-P algorithm when the EM-grid is a full Manhattan grid (see
$G_{2}$ in Figure 10), or similarly, to perform the CMEB-P algorithm when the
EM-grid is a full Euclidean grid (see $G_{3}$). In the general mixed case (see
$G_{1}$), though, all the algorithms perform decently with respect to the
optimal one.

Taking into account the Chicago map, i.e., $G_{1}$, we can immediately observe
that the adopted Matrice drone can potentially accomplish the whole mission of
up to $n=50$ deliveries with a single battery charge, regardless of the
algorithm. Specifically, on average, both OPT-P and S-OPT-P require
$29.181\mathrm{\,km}$ of travel, followed by $29.297\mathrm{\,km}$ for
CMALL-P, $29.324\mathrm{\,km}$ for CMEB-P, and finally $30.395\mathrm{\,km}$
for CEMB-P. We note that the drone, if employing Algorithm CEMB-P, may be
forced to fly beyond the aforementioned and approximated threshold of
$30\mathrm{\,km}$. However, here one can appreciate that the difference, on
average, among the best and the worst performing algorithm is $\approx
1\mathrm{\,km}$, which can be detrimental in some borderline cases. Moreover,
even neglecting the CEMB-P algorithm, the above difference is still about
$0.143\mathrm{\,km}$. Recalling that the average block length in $G_{1}$ is
$0.12\mathrm{\,km}$, this small difference could potentially result in a
delivery not being completed.

In the $G_{2}$ map in New York City, which is a full Manhattan grid, we can
safely do up to $n=35$ deliveries. In this case, we cannot rely to the CEMB-P
algorithm, since it is not suitable for full Manhattan grids. It is important
to recall that, in $G_{2}$, CMALL-P returns, by definition of the Manhattan
median, the optimal solution.

Finally, in the full Euclidean grid in Miami ($G_{3}$), due to the smaller
distances in the example, the Matrice drone can safely perform even more than
$n=100$ deliveries. Even relying to the CMEB-P algorithm, we are still able to
serve $100$ customers in the area. It is important to recall that, in $G_{3}$,
CEMB-P returns the optimal solution.

##  7 Conclusion

We considered a drone-based delivery system for the “last-mile” logistics of
small parcels, medicines, or viral tests, in EM-grids. The shortest path in an
EM-grid concatenates Euclidean and Manhattan distances. We solved the MEMP on
EM-grids whose goal is to minimize the sum of the distances between the
locations to be served and the drone’s DP. Finding the most suitable DP has
many implications that can impact the expected delivery time for customers,
the energy consumption of drones, and in general, the broader environmental
impacts, e.g., the quantity of CO2 emissions, when relying on trucks. We
propose efficient algorithms to exactly solve the problem in EM-grids in both
the full-grid and the partial-grid scenarios under the assumption that the
vertex distance is unitary. We also run our algorithms on quasi-real cases by
extracting the EM-grids from real cities in the United States.

Although we attempted to run our algorithms on real cities, we are aware that
our EM-grid model is too simplified to characterize any real-world scenario.
In future work, we intend to extend the introduced mixed-grid model to more
general layouts (e.g., a rural area inside an urban area in Central Park, New
York). Nevertheless, we believe that our work could be the starting point to
devising much more complex scenarios to better model real-world scenarios. For
example, we could use our technique of map extraction (see Section 6.3) to
construct a grid, $G$, where the each individual grid square is labeled as
“Buildings” or “Park”. A contiguous group of squares labeled “Buildings” is
effectively a Manhattan grid, and a contiguous group of squares labeled “Park”
is a Euclidean grid. Then, we could sub-divide $G$ into multiple EM-grids, and
apply our algorithms suitably in each EM-grid.

Another interesting variant to study is the use of multiple drones that are
responsible for delivering packages to different partitions of the customers.
In this case, we will need efficient clustering algorithms for minimizing
suitable metrics.

## References

  * [1] J. Wivou, L. Udawatta, A. Alshehhi, E. Alzaabi, A. Albeloshi, and S. Alfalasi, “Air quality monitoring for sustainable systems via drone based tech.” in _Intl. Conf. on Information and Automation for Sustainability (ICIAfS)_. IEEE, 2016, pp. 1–5.
  * [2] T. Calamoneri, F. Corò, and S. Mancini, “A realistic model to support rescue operations after an earthquake via uavs,” _IEEE Access_ , 2022.
  * [3] C. Y. Justin, A. P. Payan, and D. N. Mavris, “Integrated fleet assignment and scheduling for environmentally friendly electrified regional air mobility,” _Transportation Research Part C: Emerging Technologies_ , vol. 138, p. 103567, 2022.
  * [4] C. Qu, F. Betti Sorbelli, R. Singh, P. Calyam, and S. K. Das, “Environmentally-aware and energy-efficient multi-drone coordination and networking for disaster response,” _IEEE Transactions on Network and Service Management_ , 2023.
  * [5] X. Li, D. Guo, H. Yin, and G. Wei, “Drone-assisted public safety wireless broadband network,” in _2015 IEEE Wireless Communications and Networking Conf._ IEEE, 2015, pp. 323–328.
  * [6] D. He, S. Chan, and M. Guizani, “Drone-assisted public safety networks: The security aspect,” _IEEE Communications Magazine_ , vol. 55, no. 8, pp. 218–223, 2017.
  * [7] Y. V. Pant, M. Z. Li, A. Rodionova, R. A. Quaye, H. Abbas, M. S. Ryerson, and R. Mangharam, “Fads: A framework for autonomous drone safety using temporal logic-based trajectory planning,” _Transportation Research Part C: Emerging Technologies_ , vol. 130, p. 103275, 2021.
  * [8] V. Niculescu, D. Palossi, M. Magno, and L. Benini, “Energy-efficient, precise uwb-based 3-d localization of sensor nodes with a nano-uav,” _IEEE Internet of Things Journal_ , 2022.
  * [9] T. Moribe, H. Okada, K. Kobayashl, and M. Katayama, “Combination of a wireless sensor network and drone using infrared thermometers for smart agriculture,” in _2018 15th IEEE Annual Consumer Communications & Networking Conf. (CCNC)_. IEEE, 2018, pp. 1–2.
  * [10] A. M. Jawad, H. M. Jawad, R. Nordin, S. K. Gharghan, N. F. Abdullah, and M. J. Abu-Alshaeer, “Wireless power transfer with magnetic resonator coupling and sleep/active strategy for a drone charging station in smart agriculture,” _IEEE Access_ , vol. 7, pp. 139 839–139 851, 2019.
  * [11] H. P. D. Nguyen and D. D. Nguyen, “Drone application in smart cities: The general overview of security vulnerabilities and countermeasures for data communication,” _Development and Future of Internet of Drones_ , pp. 185–210, 2021.
  * [12] T. Preethika, P. Vaishnavi, J. Agnishwar, K. Padmanathan, S. Umashankar, S. Annapoorani, M. Subash, and K. Aruloli, “Artificial intelligence and drones to combat covid-19,” _Preprints_ , 2020.
  * [13] D. S. Jat and C. Singh, “Artificial intelligence-enabled robotic drones for covid-19 outbreak,” in _Intelligent Systems and Methods to Combat Covid-19_. Springer, 2020, pp. 37–46.
  * [14] S. Anggraeni, A. Maulidina, M. W. Dewi, S. Rahmadianti, Y. P. C. Rizky, Z. F. Arinalhaq, D. Usdiyana, A. B. D. Nandiyanto, A. S. Al-Obaidi _et al._ , “The deployment of drones in sending drugs and patient blood samples covid-19,” _Indonesian Journal of Science and Technology_ , pp. 18–25, 2020\.
  * [15] V. Pandit and A. Poojari, “A study on amazon prime air for feasibility and profitability: A graphical data analysis,” _IOSR Journal of Business and Management_ , vol. 16, no. 11, pp. 06–11, 2014.
  * [16] A. Welch, “A cost-benefit analysis of amazon prime air,” 2015.
  * [17] S. M. Shavarani, M. G. Nejad, F. Rismanchian, and G. Izbirak, “Application of hierarchical facility location problem for optimization of a drone delivery system: a case study of amazon prime air in the city of san francisco,” _The International Journal of Advanced Manufacturing Technology_ , vol. 95, pp. 3141–3153, 2018.
  * [18] J. Pepitone, “Domino’s tests drone pizza delivery,” _CNNMoney, June_ , vol. 4, 2013.
  * [19] Q. M. Ha, Y. Deville, Q. D. Pham, and M. H. Hà, “Heuristic methods for the traveling salesman problem with drone,” _Comput. Sci_ , 2015.
  * [20] S. Sawadsitang, D. Niyato, P. S. Tan, P. Wang, and S. Nutanong, “Multi-objective optimization for drone delivery,” in _90th Vehicular Tech. Conf. (VTC2019-Fall)_. IEEE, 2019, pp. 1–5.
  * [21] I. Dayarian, M. Savelsbergh, and J.-P. Clarke, “Same-day delivery with drone resupply,” _Transp. Science_ , vol. 54, no. 1, pp. 229–249, 2020.
  * [22] N. Agatz, P. Bouman, and M. Schmidt, “Optimization approaches for the traveling salesman problem with drone,” _Transp. Science_ , 2018.
  * [23] FAA, “Faadronezone,” https://faadronezone.faa.gov, 2023, (Accessed on May 2023).
  * [24] M. Kunovjanek and C. Wankmüller, “Containing the covid-19 pandemic with drones-feasibility of a drone enabled back-up transport system,” _Transport Policy_ , vol. 106, pp. 141–152, 2021.
  * [25] L. Bartoli, F. Betti Sorbelli, F. Corò, C. M. Pinotti, and A. Shende, “Exact and approximate drone warehouse for a mixed landscape delivery system,” in _Intl. Conf. on Smart Computing (SMARTCOMP)_. IEEE, 2019, pp. 266–273.
  * [26] F. Betti Sorbelli, F. Corò, C. M. Pinotti, and A. Shende, “Automated picking system employing a drone,” in _Intl. Conf. on Distributed Computing in Sensor Systems (DCOSS)_. IEEE, 2019, pp. 633–640.
  * [27] C. C. Murray and A. G. Chu, “The flying sidekick traveling salesman problem: Optimization of drone-assisted parcel delivery,” _Transp. Research Part C: Emerging Tech._ , vol. 54, pp. 86–109, 2015.
  * [28] C. C. Murray and R. Raj, “The multiple flying sidekicks traveling salesman problem: Parcel delivery with multiple drones,” _Transp. Research Part C: Emerging Tech._ , vol. 110, pp. 368–398, 2020.
  * [29] M. Dell’Amico, R. Montemanni, and S. Novellani, “Exact models for the flying sidekick traveling salesman problem,” _Intl. Trans. in Operational Research_ , vol. 29, no. 3, pp. 1360–1393, 2022.
  * [30] K. Kloster, M. Moeini, D. Vigo, and O. Wendt, “The multiple traveling salesman problem in presence of drone-and robot-supported packet stations,” _European Journal of Operational Research_ , vol. 305, no. 2, pp. 630–643, 2023.
  * [31] D. Schermer, M. Moeini, and O. Wendt, “A hybrid vns/tabu search algorithm for solving the vehicle routing problem with drones and en route operations,” _Computers & Operations Research_, vol. 109, pp. 134–158, 2019.
  * [32] J. R. Brown, M. A. Bushuev, and A. L. Guiffrida, “Distance metrics matter: analysing optimisation algorithms for the last mile problem,” _Intl. Journal of Logistics Systems and Management_ , vol. 38, no. 2, pp. 151–174, 2021\.
  * [33] N. Ai, J. Zheng, X. Chen, and K. Kawamura, “Neighborhood-specific traffic impact analysis of restaurant meal delivery trips: Planning implications and case studies in chicago,” _Journal of Urban Planning and Development_ , vol. 147, no. 2, p. 05021013, 2021.
  * [34] S. Poikonen, B. Golden, and E. A. Wasil, “A branch-and-bound approach to the traveling salesman problem with a drone,” _INFORMS Journal on Computing_ , vol. 31, no. 2, pp. 335–346, 2019.
  * [35] F. Betti Sorbelli, F. Corò, S. K. Das, L. Palazzetti, and C. M. Pinotti, “On the scheduling of conflictual deliveries in a last-mile delivery scenario with truck-carried drones,” _Pervasive and Mobile Computing_ , vol. 87, p. 101700, 2022.
  * [36] F. Betti Sorbelli, F. Corò, L. Palazzetti, C. M. Pinotti, and G. Rigoni, “How the wind can be leveraged for saving energy in a truck-drone delivery system,” _IEEE Transactions on Intelligent Transportation Systems_ , 2023\.
  * [37] Y. Li, W. Yang, and B. Huang, “Impact of uav delivery on sustainability and costs under traffic restrictions,” _Mathematical Problems in Engineering_ , vol. 2020, 2020.
  * [38] M. Salama and S. Srinivas, “Joint optimization of customer location clustering and drone-based routing for last-mile deliveries,” _Transp. Research Part C: Emerging Tech._ , vol. 114, pp. 620–642, 2020.
  * [39] O. Dukkanci, B. Y. Kara, and T. Bektaş, “Minimizing energy and cost in range-limited drone deliveries with speed optimization,” _Transp. Research Part C: Emerging Tech._ , vol. 125, p. 102985, 2021.
  * [40] R. Pinto and A. Lagorio, “Point-to-point drone-based delivery network design with intermediate charging stations,” _Transp. Research Part C: Emerging Tech._ , vol. 135, p. 103506, 2022.
  * [41] A. Karak and K. Abdelghany, “The hybrid vehicle-drone routing problem for pick-up and delivery services,” _Transp. Research Part C: Emerging Tech._ , vol. 102, pp. 427–449, 2019.
  * [42] G. Clarke and J. W. Wright, “Scheduling of vehicles from a central depot to a number of delivery points,” _Operations research_ , vol. 12, no. 4, pp. 568–581, 1964.
  * [43] I. Hong, M. Kuby, and A. T. Murray, “A range-restricted recharging station coverage model for drone delivery service planning,” _Transp. Research Part C: Emerging Tech._ , vol. 90, pp. 198–212, 2018.
  * [44] T. Yamaguchi and I. Kaji, “Some location problems on grid graphs,” _Electronics and Communications in Japan (Part I: Communications)_ , vol. 70, no. 10, pp. 31–40, 1987.
  * [45] Y. Vardi and C.-H. Zhang, “The multivariate l1-median and associated data depth,” _Proceedings of the National Academy of Sciences_ , vol. 97, no. 4, pp. 1423–1426, 2000.
  * [46] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, _Introduction to algorithms_. MIT press, 2009.
  * [47] DJI, _Matrice 300 RTK - Specs_ , 2022, https://www.dji.com/matrice-300/specs.