The Drone Delivery Problem#

Optimization Model for Drone Delivery

This notebooks outlines the optimization model designed to determine the minimum distance for drone delivery routing problems. The model employs a Mixed Integer Linear Programming (MILP) approach to minimize the total travel distance required to deliver packages from a warehouse to various delivery points.

Libraries

To begin, we import the necessary libraries.

pip install quantagonia
pip install folium
API_KEY="Your-API-KEY" # if you don't have one, head over to https://platform.quantagonia.com/ (free account available)
import numpy as np
import pandas as pd
import time as ttt
import random
import math
import folium

from pulp import *
from quantagonia import HybridSolverParameters
import quantagonia.mip.pulp_adapter as pulp_adapter

Problem Parameters and Data Preparation

Next, we set up the initial parameters.

# Define problem parameters
NUM_LOCATIONS = 12 # Number of Location Points Including Warehouse
DRONE_CAPACITY = 3

Generate coordinates for the selected area, which in our case is Munich. To gain an initial perspective on our problem, we will visualize the data using the Folium.

random.seed(2)

# Latitude and Longitude ranges for Munich
lat_range = (48.114972911995785, 48.17622969284977) # Approximate latitude range for Munich
lon_range = (11.536111191947134, 11.611750294370422) # Approximate longitude range for Munich

# Generate random coordinates within the specified range
coordinates=np.zeros(NUM_LOCATIONS,dtype=object)
for i in range(NUM_LOCATIONS):
    coordinates[i] = (round(random.uniform(lat_range[0], lat_range[1]), 17), round(random.uniform(lon_range[0], lon_range[1]), 8))

    # Create a map centered around the average coordinates
avg_lat = sum(lat for lat, lon in coordinates) / len(coordinates)
avg_lon = sum(lon for lat, lon in coordinates) / len(coordinates)
map_munich = folium.Map(location=[avg_lat, avg_lon], zoom_start=10)

# Add markers for each coordinate
for i, coord in enumerate(coordinates):
    if i == 0:
        folium.Marker(location=coord, popup=f"Name: Warehouse, Lat: {coord[0]}, Lon: {coord[1]}", icon=folium.Icon(color='green')).add_to(map_munich)
    else:
        folium.Marker(location=coord, popup=f"Name: Delivery Point {i}, Lat: {coord[0]}, Lon: {coord[1]}").add_to(map_munich)

# Save the map to an HTML file
map_munich.save("Munich_coordinates_map.html")

# If running in a Jupyter notebook, display the map
map_munich

Distance Calculation

We use the Haversine formula to calculate the distance between two geographic coordinates and then construct the distance matrix for the model.

# Haversine Distance Function
def haversine_distance(coord1, coord2):
    R = 6371  # Radius of the Earth in kilometers
    lat1, lon1 = math.radians(coord1[0]), math.radians(coord1[1])
    lat2, lon2 = math.radians(coord2[0]), math.radians(coord2[1])
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    a = math.sin(dlat / 2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2)**2
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
    return R * c

# Construct the distance matrix
distance_matrix = np.zeros((len(coordinates),len(coordinates)))
for i in range(len(coordinates)):
    for j in range(len(coordinates)):
        distance_matrix[i,j] = haversine_distance(coordinates[i], coordinates[j])

Sets and Indices:

Now, we define the sets:

  • \(i,j \in \text{locationset} = \{0,1,\ldots,(n-1)\}\), the set of locations (including warehouse, \(idx=0\))

  • \(i,j \in \text{locationset_without_wh} = \{1,2,\ldots,(n-1)\}\), the set of locations (excluding warehouse, \(idx=0\))

  • \(t \in \text{tours} = \{0,1,\ldots,\text{Maximum_Tours}\}\), the set of tours

# SETS
locationset = list(range(NUM_LOCATIONS)) ## Location set includinng warehouse (idx=0)
locationset_without_wh = list(range(1,NUM_LOCATIONS)) ## Location set excluding warehouse (idx=0)
tours = list(range(math.ceil((NUM_LOCATIONS-1)/DRONE_CAPACITY))) ## Tour set

Decision Variables

We define the decision variables for the model to be:

  • \(x_{ijt}\): binary Variable, 1 if location \(j\) is visited after location \(i\) within tour \(t\), otherwise 0.

  • \(y_{it}\): binary Variable, 1 if location \(i\) is visited in tour \(t\), otherwise \(0\).

  • \(u_i\): auxiliary integer variable to eliminate subtours, representing the order in which delivery point \(i\) is visited.

#DECISION VARIABLES

### x[i][j][t] - Binary: if location j is visited after location i within tour t, then 1; otherwise 0
x = []
for i in locationset:
    x.append([])
    for j in locationset:
        x[i].append([])
        for t in tours:
            x[i][j].append(LpVariable(name = 'x(%d,%d,%d)' % (i,j,t), cat='Binary'))
### y[i][t] - Binary: if location i is visited in tour t, then 1; otherwise 0
y = []
for i in locationset:
    y.append([])
    for t in tours:
        y[i].append(LpVariable(name = 'y(%d,%d)' % (i,t), cat='Binary'))

### u[i] - Auxiliary variable to eliminate subtours, representing the order in which delivery point i is visited
u = []
for i in locationset:
    u.append(LpVariable(name = 'u(%d)' % (i),lowBound=0, cat='Integer'))

Objective Function

We define the objective function to minimize the total travel distance:

\begin{equation} \text{Minimize} \quad Z = \sum_{i \in \text{locationset}} \sum_{j \in \text{locationset}} \sum_{t \in \text{tours}} x_{ijt} \cdot \text{distance_matrix}_{ij} \tag{0} \end{equation}

Z = LpProblem("Deliver_Packages", LpMinimize)
### OBJECTIVE FUNCTION
Z += lpSum(x[i][j][t]*distance_matrix[i][j] for i in locationset for j in locationset for t in tours)

Constraints

Constraint 1:

The number of entries and exits for each location must be the same.

\[\begin{split}\begin{align} \sum_{i \in \text{locationset}} x_{ijt} &= \sum_{i \in \text{locationset}} x_{jit}, \quad \forall j \in \text{locationset}, \forall t \in \text{tours} \tag{1}\\ \end{align}\end{split}\]

Constraint 2a:

The number of entries to each location \(j\) must be equal to 1 if it is visited in tour \(t\).

\[\begin{split}\begin{align} \sum_{i \in \text{locationset}, i \neq j} x_{ijt} &= y_{jt}, \quad \forall j \in \text{locationset}, \forall t \in \text{tours} \tag{2a}\\ \end{align}\end{split}\]

Constraint 2b:

The number of departures from each location \(i\) must be equal to 1 if it is visited in tour \(t\).

\[\begin{split}\begin{align} \sum_{j \in \text{locationset}, j \neq i} x_{ijt} &= y_{it}, \quad \forall i \in \text{locationset}, \forall t \in \text{tours} \tag{2b}\\ \end{align}\end{split}\]

Constraint 3:

Each delivery location must be visited.

\[\begin{split}\begin{align} \sum_{t \in \text{tours}} y_{it} &= 1, \quad \forall i \in \text{locationset_without_wh} \tag{3}\\ \end{align}\end{split}\]

Constraint 4:

Capacity constraint: The maximum number of customers visited is limited by the drone’s capacity. Additionally, the warehouse node must be visited if any tour is activated.

\[\begin{split}\begin{align} \sum_{i \in \text{locationset_without_wh}} y_{it} &\leq \text{DRONE_CAPACITY}*y_{0t}, \quad \forall t \in \text{tours} \tag{4}\\ \end{align}\end{split}\]

Constraint 5a:

Subtour elimination and order constraints (Miller-Tucker-Zemlin (MTZ) formulation).

\[\begin{split}\begin{align} u_i - u_j + \text{DRONE_CAPACITY} \cdot x_{ijt} &\leq \text{DRONE_CAPACITY} - 1, \quad \forall i, j \in \text{locationset_without_wh}, i \neq j, \forall t \in \text{tours} \tag{6a}\\ \end{align}\end{split}\]

Constraint 5b:

Subtour elimination and order constraints (MTZ formulation).

\[\begin{split}\begin{align} 1 \leq u_i &\leq \text{DRONE_CAPACITY}, \quad \forall i \in \text{locationset_without_wh} \tag{6b}\\ \end{align}\end{split}\]
### CONSTRAINTS
# Constraint 1: The number of entries and exits for each location must be the same.
for j in locationset:
    for t in tours:
        Z += lpSum(x[i][j][t] for i in locationset) == lpSum(x[j][i][t] for i in locationset)

# Constraint 2a: The number of entries to each location j must be equal to 1 if it is visited in tour t.
for j in locationset:
    for t in tours:
        Z +=lpSum(x[i][j][t] for i in locationset if i!=j) == y[j][t]

# Constraint 2b: The number of departures from each location i must be equal to 1 if it is visited in tour t.
for i in locationset:
    for t in tours:
        Z +=lpSum(x[i][j][t] for j in locationset if j!=i) == y[i][t]

# Constraint 3: Each delivery location must be visited.
for i in locationset_without_wh:
    Z += lpSum(y[i][t] for t in tours) == 1

# Constraint 4: Capacity constraint: The maximum number of customers visited is limited by the drone's capacity.
# Additionally, the warehouse node must be visited if any tour is activated.
for t in tours:
    Z += lpSum(y[i][t] for i in locationset_without_wh) <= DRONE_CAPACITY*y[0][t]

# Constraint 5a: Subtour Elimination and Order Constraints (Miller-Tucker-Zemlin (MTZ) formulation).
for t in tours:
    for i in locationset_without_wh:
        for j in locationset_without_wh:
            if i==j:
                continue
            Z += (u[i]-u[j]+DRONE_CAPACITY*x[i][j][t] <= DRONE_CAPACITY-1)

# Constraint 5b: Subtour Elimination and Order Constraints (MTZ formulation).
for t in tours:
    for i in locationset_without_wh:
        Z += (u[i] <= DRONE_CAPACITY)
        Z += (u[i] >= 1)

Solving the Model

We define the parameters (e.g. Time, Gap) and solve the model using a the HybridSolver from Quantagonia.

params = HybridSolverParameters()
params.set_time_limit(600)
params.set_seed(0)
q_solver = pulp_adapter.HybridSolver_CMD(api_key=API_KEY, params=params)
status = Z.solve(solver=q_solver)

print("Status:",LpStatus[status],"\nObjective Function Value: %.2f"%Z.objective.value(),"km")
✔ Queued job with jobid 0028b7fb-9985-4f49-a3dc-f18e118b3b0a...
✔ Job 0028b7fb-9985-4f49-a3dc-f18e118b3b0a unqueued, processing...

Quantagonia HybridSolver version 1.1.1781
Copyright (c) 2024 Quantagonia GmbH.
HybridSolver integrates various open-source packages; see release notes.

User-specified parameters:
Set parameter 'time_limit' to value '600.0'.
Set parameter 'seed' to value '0'.

Read b145cc629f4d4258ae4d044a766029be-pulp.mps in 0.59s.
Minimize a MILP with 687 constraints and 635 variables (624 binary, 11 integer, 0 implied integer,0 continuous).

Presolving model. Presolved model in 0.0s.
Reduced model has 551 constraints and 576 variables (565 binary, 11 integer, 0 continuous).


------------------------------------------------------------------------
       Nodes |      Incumbent |          Bound |   Gap (%) |  Time (s) |
------------------------------------------------------------------------
           1 |            inf |     0.00000000 |       inf |      0.00 |
 *         1 |     46.7744772 |     19.8227293 |     57.62 |      0.00 |
 *         3 |     43.9381887 |     27.2419397 |     38.00 |      0.02 |
 *        59 |     43.3015915 |     27.2488165 |     37.07 |      0.43 |
 *        82 |     43.0262083 |     27.2488165 |     36.67 |      0.68 |
 *       162 |     42.3020208 |     27.2496758 |     35.58 |      1.27 |
 *       229 |     42.1412607 |     27.2605539 |     35.31 |      1.80 |
         386 |     42.1412607 |     27.2622454 |     35.31 |      2.93 |
 *       539 |     42.1061606 |     27.2622454 |     35.25 |      4.23 |
 *       880 |     41.9052291 |     28.0060363 |     33.17 |      4.76 |
 *       884 |     41.8789535 |     28.0060363 |     33.13 |      4.86 |
 *      1131 |     41.7559333 |     29.8798681 |     28.44 |      5.69 |
        2017 |     41.7559333 |     31.0858794 |     25.55 |     11.48 |
        2142 |     41.7559333 |     31.4118209 |     24.77 |     12.80 |
        2338 |     41.7559333 |     31.6075000 |     24.30 |     14.04 |
        2521 |     41.7559333 |     31.6556399 |     24.19 |     15.06 |
        2770 |     41.7559333 |     32.1571995 |     22.99 |     16.10 |
        2878 |     41.7559333 |     32.4013208 |     22.40 |     17.15 |
        3036 |     41.7559333 |     32.6622455 |     21.78 |     18.62 |
        3249 |     41.7559333 |     32.7038099 |     21.68 |     19.90 |
        3465 |     41.7559333 |     32.8898708 |     21.23 |     21.54 |
        3682 |     41.7559333 |     32.9232617 |     21.15 |     22.86 |
        3908 |     41.7559333 |     33.0918351 |     20.75 |     24.07 |
        4073 |     41.7559333 |     33.2434330 |     20.39 |     25.97 |
        4385 |     41.7559333 |     33.3867754 |     20.04 |     27.48 |
------------------------------------------------------------------------
       Nodes |      Incumbent |          Bound |   Gap (%) |  Time (s) |
------------------------------------------------------------------------
        5280 |     41.7559333 |     33.5461010 |     19.66 |     28.54 |
        5425 |     41.7559333 |     33.7102887 |     19.27 |     29.89 |
        5811 |     41.7559333 |     33.7967712 |     19.06 |     31.01 |
        6093 |     41.7559333 |     33.8393989 |     18.96 |     32.40 |
        6491 |     41.7559333 |     34.0827942 |     18.38 |     33.42 |
        6590 |     41.7559333 |     34.2933256 |     17.87 |     34.68 |
        6883 |     41.7559333 |     34.4214993 |     17.57 |     35.74 |
        7370 |     41.7559333 |     34.6465789 |     17.03 |     37.63 |
        7620 |     41.7559333 |     34.8405817 |     16.56 |     38.85 |
        8069 |     41.7559333 |     34.9895971 |     16.20 |     40.07 |
        8380 |     41.7559333 |     35.1596993 |     15.80 |     41.13 |
        8557 |     41.7559333 |     35.5214184 |     14.93 |     42.23 |
        8926 |     41.7559333 |     35.6900359 |     14.53 |     43.50 |
        9366 |     41.7559333 |     35.8195714 |     14.22 |     44.66 |
        9814 |     41.7559333 |     35.8797123 |     14.07 |     45.72 |
       10041 |     41.7559333 |     36.0815336 |     13.59 |     46.83 |
       10391 |     41.7559333 |     36.3871303 |     12.86 |     48.00 |
       10626 |     41.7559333 |     36.5415914 |     12.49 |     49.02 |
       10924 |     41.7559333 |     36.7241490 |     12.05 |     50.25 |
       11593 |     41.7559333 |     36.9130140 |     11.60 |     51.30 |
       11826 |     41.7559333 |     37.0363145 |     11.30 |     52.47 |
       12032 |     41.7559333 |     37.1233716 |     11.09 |     53.69 |
       12624 |     41.7559333 |     37.4367840 |     10.34 |     55.43 |
       13483 |     41.7559333 |     37.6922509 |      9.73 |     56.55 |
       13823 |     41.7559333 |     37.9103551 |      9.21 |     57.75 |
------------------------------------------------------------------------
       Nodes |      Incumbent |          Bound |   Gap (%) |  Time (s) |
------------------------------------------------------------------------
       14282 |     41.7559333 |     38.2179958 |      8.47 |     58.94 |
       14597 |     41.7559333 |     38.4545453 |      7.91 |     59.95 |
       14974 |     41.7559333 |     38.6770731 |      7.37 |     60.99 |
       15903 |     41.7559333 |     38.9883507 |      6.63 |     62.13 |
       16327 |     41.7559333 |     39.2174342 |      6.08 |     63.20 |
       16788 |     41.7559333 |     39.4652921 |      5.49 |     64.27 |
       17336 |     41.7559333 |     39.7677276 |      4.76 |     65.38 |
       18056 |     41.7559333 |     40.2870297 |      3.52 |     66.44 |
       18918 |     41.7559333 |     41.7534278 |      0.01 |     67.44 |
------------------------------------------------------------------------

Optimal solution found (within relative tolerance 0.01%).

Solver Results:
 - Solution Status: Optimal
 - Wall Time: 67.4439 seconds
 - Objective: 41.7559333
 - Bound: 41.7534278
 - Absolute Gap: 0.0025
 - Relative Gap: 0.006%
 - Nodes: 18918
 - Best solution found at node 1131 after 6.2713 seconds
Finished processing job 0028b7fb-9985-4f49-a3dc-f18e118b3b0a...
Optimal
Status: Optimal
Objective Function Value: 41.76 km

Results

We extract the results:

# GET THE END RESULTS OF THE VARIABLES

resultx = np.zeros((len(locationset),len(locationset),len(tours)))
for i in locationset:
    for j in locationset:
        for t in tours:
            resultx[i,j,t] = x[i][j][t].varValue

resulty = np.zeros((len(locationset),len(tours)))
for i in locationset:
    for t in tours:
         resulty[i,t] = y[i][t].varValue

## Extract the sequences
result_u = np.zeros(len(locationset))
for i in locationset_without_wh:
    result_u[i]=u[i].varValue

## Extract the squences of the locations and which in tour they are visited
sequence_tour = np.zeros((NUM_LOCATIONS,2))
sequence_tour[:,1] = -1
for i in locationset_without_wh:
    sequence_tour[i,0] = int(round(result_u[i]))
    for t in tours:
        if int(round(resulty[i,t]))==1:
            sequence_tour[i,1] = t

# Construct the Output DataFrame
output_df = pd.DataFrame(sequence_tour[1:], columns=['Sequence','Tour'])
output_df['Location_Index'] = range(1,len(sequence_tour))
output_df['Coordinates'] = coordinates[1:]
output_df.sort_values(['Tour', 'Sequence'], ascending=[True, True],inplace=True)
for i in range(1,len(output_df)):
    if output_df.iloc[i,1] == output_df.iloc[i-1,1]:
        output_df.iloc[i,0] = output_df.iloc[i-1,0] + 1
    else:
        output_df.iloc[i,0] = 1

Visualizing the Results

We visualize the results using Folium:

### MAP VIEW OF THE SOLUTION

# Define the list of colors for the tours (currently max tours is limited by 15)
colors = ['purple', 'red', 'blue', 'black', 'darkred', 'lightred', 'beige',
          'darkblue', 'cadetblue', 'darkpurple', 'pink', 'lightblue',
          'gray', 'orange', 'lightgray']

# Calculate the average coordinates for centering the map
avg_lat = sum(lat for lat, lon in coordinates) / len(coordinates)
avg_lon = sum(lon for lat, lon in coordinates) / len(coordinates)

# Create a map centered around the average coordinates
map_munich = folium.Map(location=[avg_lat, avg_lon], zoom_start=10)

# Add the first marker as the warehouse with a green icon
folium.Marker(
    location=coordinates[0],
    popup=f"Name: Warehouse, Lat: {coordinates[0][0]}, Lon: {coordinates[0][1]}",
    icon=folium.Icon(color='green')).add_to(map_munich)

# Add markers for each coordinate
for i in range(len(output_df)):

    # Add markers for other locations with DivIcon displaying the sequence number
    folium.Marker(
        location=output_df['Coordinates'].iloc[i],
        popup=f"Location {output_df['Location_Index'].iloc[i]}: Lat: {output_df['Coordinates'].iloc[i][0]}, Lon: {output_df['Coordinates'].iloc[i][1]}",
        icon=folium.DivIcon(html=f'<div style="font-size: 16pt; color : black">{int(output_df["Sequence"].iloc[i])}</div>')
    ).add_to(map_munich)

    # Add circle markers for other locations with specified colors
    folium.CircleMarker(
        location=output_df['Coordinates'].iloc[i],
        radius=5,
        color=colors[int(output_df['Tour'].iloc[i])%len(colors)],
        fill=True,
        fill_color=colors[int(output_df['Tour'].iloc[i])%len(colors)]
    ).add_to(map_munich)

# Create a dictionary to store coordinates by tour
tour_coordinates = {}
for i in range(len(output_df)):
    tour = output_df['Tour'].iloc[i]
    if tour not in tour_coordinates:
        tour_coordinates[tour] = []
    tour_coordinates[tour].append(output_df['Coordinates'].iloc[i])

# Draw dashed lines for each tour
for tour, coords in tour_coordinates.items():
    # Connect warehouse to the first location in the tour
    folium.PolyLine(
        locations=[coordinates[0], coords[0]],
        color=colors[int(tour)%len(colors)],
        weight=2,
        opacity=0.5,
        dash_array='1, 1'
    ).add_to(map_munich)

    # Connect each location in the tour by its sequence
    folium.PolyLine(
        locations=coords,
        color=colors[int(tour)%len(colors)],
        weight=2,
        opacity=0.5,
        dash_array='1, 1'
    ).add_to(map_munich)

    # Connect the last location in a tour back to the warehouse
    folium.PolyLine(
        locations=[coords[-1], coordinates[0]],
        color=colors[int(tour)%len(colors)],
        weight=2,
        opacity=0.5,
        dash_array='1, 1'
    ).add_to(map_munich)

map_munich.save("map_munich.html")
map_munich

Need help with modeling? We are happy to coach you through your model formulation. Reach out to us at help@quantagonia.com or https://www.quantagonia.com/contact.