Classic VRPs¶

This notebook shows how to use PyVRP to solve two classic variants of the VRP: the capacitated vehicle routing problem (CVRP), and the vehicle routing problem with time windows (VRPTW). It builds on the tutorial by solving much larger instances, and going into more detail about the various plotting tools and diagnostics available in PyVRP.

A CVRP instance is defined on a complete graph \(G=(V,A)\), where \(V\) is the vertex set and \(A\) is the arc set. The vertex set \(V\) is partitioned into \(V=\{0\} \cup V_c\), where \(0\) represents the depot and \(V_c=\{1, \dots, n\}\) denotes the set of \(n\) customers. Each arc \((i, j) \in A\) has a weight \(d_{ij} \ge 0\) that represents the travel distance from \(i \in V\) to \(j \in V\). Each customer \(i \in V_c\) has a demand \(q_{i} \ge 0\). The objective is to find a feasible solution that minimises the total distance.

A VRPTW instance additionally incorporates time aspects into the problem. For the sake of exposition we assume the travel duration \(t_{ij} \ge 0\) is equal to the travel distance \(d_{ij}\) in this notebook. Each customer \(i \in V_c\) has a service time \(s_{i} \ge 0\) and a (hard) time window \(\left[e_i, l_i\right]\) that denotes the earliest and latest time that service can start. A vehicle is allowed to arrive at a customer location before the beginning of the time window, but it must wait for the window to open to start the delivery. Each vehicle must return to the depot before the end of the depot time window \(H\). The objective is to find a feasible solution that minimises the total distance.

Let’s first import what we will use in this notebook.

[1]:

import matplotlib.pyplot as plt
from tabulate import tabulate
from vrplib import read_solution

from pyvrp import read, solve
from pyvrp.plotting import (
    plot_coordinates,
    plot_instance,
    plot_result,
    plot_route_schedule,
)
from pyvrp.stop import MaxIterations, MaxRuntime

The capacitated VRP¶

Reading the instance¶

We will solve the X-n439-k37 instance, which is part of the X instance set that is widely used to benchmark CVRP algorithms. The function pyvrp.read reads the instance file and converts it to a ProblemData instance. We pass the argument round_func="round" to compute the Euclidean distances rounded to the nearest integral, which is the convention for the X benchmark set. We also load the best known solution to evaluate our solver later on.

[2]:

INSTANCE = read("data/X-n439-k37.vrp", round_func="round")
BKS = read_solution("data/X-n439-k37.sol")

Let’s plot the instance and see what we have.

[3]:

_, ax = plt.subplots(figsize=(8, 8))
plot_coordinates(INSTANCE, ax=ax)
plt.tight_layout()

Solving the instance¶

We will use the default solve method that PyVRP provides for data instances.

[4]:

result = solve(INSTANCE, stop=MaxIterations(2000), seed=42, display=False)
print(result)

Solution results
================
    # routes: 37
     # trips: 37
   # clients: 438
   objective: 36737
    distance: 36737
    duration: 36737
# iterations: 2000
    run-time: 44.70 seconds

Routes
------
Route #1: 414 217 236 59 105 434 311 8 2 169 3 348
Route #2: 411 410 97 400 267 386 299 223 349 425 133 370
Route #3: 326 26 260 92 275 41 406 308 202 149 172
Route #4: 346 211 347 206 57 296 375 43 281 239 42 335
Route #5: 155 71 228 162 435 166 345 385 438 381 404 270
Route #6: 72 218 392 139 200 145 122 366 384 403 428 402
Route #7: 249 280 312 250 421 409 110 303 245 241 360 227
Route #8: 380 377 433 337 242 342 221 393 237 325 121 195
Route #9: 418 416 407 412 17 89 293 339 66 83 388 225
Route #10: 437 233 324 229 268 115 353 422 44
Route #11: 329 319 309 297 126 86 315 352 264 372 391 423
Route #12: 243 321 413 420 396 432 351 266 289 253 383 323
Route #13: 257 271 338 252 101 285 193 215 153 15 159 65
Route #14: 283 246 138 376 47 286 350 341 98 137 251 146
Route #15: 204 176 131 6 56 22 210 88 173 118 91 154
Route #16: 79 189 61 130 80 344 7 197 431 144 334 25
Route #17: 62 371 109 90 113 116 274 196 75 140 1 367
Route #18: 207 397 390 184 395 287 387 134 426 5 31 333
Route #19: 430 248 19 330 108 340 240 220 331 361 161 58
Route #20: 28 343 389 327 255 302 408 152 177 135 157 24
Route #21: 132 4 21 378 112 16 67 230 34 73 175 401
Route #22: 244 171 13 103 117 63 74 354 209 11 181 84
Route #23: 265 364 259 30 291 160 292 174 382 192 235 165
Route #24: 313 178 310 222 125 190 356 290 368 300 424 417
Route #25: 363 357 214 198 272 123 64 10 37 301 379 394
Route #26: 81 179 96 256 78 182 216 106 399 77 99 124
Route #27: 284 320 234 405 168 95 52 328 213 279 114 316
Route #28: 170 183 142 219 205 150 282 224 273 269 120 332
Route #29: 39 232 35 141 208 374 322 304 419 336 369 188
Route #30: 53 163 151 199 54 27 100 51 164 186 128 20
Route #31: 107 306 314 148 68 69 85 12 48 180
Route #32: 38 167 203 247 40 46 212 94 49 45 277 111
Route #33: 201 76 119 127 102 87 29 158 191 129 261 262
Route #34: 194 50 18 307 14 436 373 317 398 36 427 156
Route #35: 33 23 318 294 185 143 147 60 104 70 258 93
Route #36: 254 9 429 355 263 55 288 298 278 226 276 415
Route #37: 358 295 359 365 238 305 82 32 362 187 136 231

[5]:

gap = 100 * (result.cost() - BKS["cost"]) / BKS["cost"]
print(f"Found a solution with cost: {result.cost()}.")
print(f"This is {gap:.1f}% worse than the best known", end=" ")
print(f"solution, which is {BKS['cost']}.")

Found a solution with cost: 36737.
This is 1.0% worse than the best known solution, which is 36391.

We’ve managed to find a very good solution quickly!

The Result object also contains useful statistics about the optimisation. We can now plot these statistics as well as the final solution use plot_result.

[6]:

fig = plt.figure(figsize=(15, 9))
plot_result(result, INSTANCE, fig)
fig.tight_layout()

PyVRP internally uses a genetic algorithm consisting of a population of feasible and infeasible solutions. These solutions are iteratively combined into new offspring solutions, that should result in increasingly better solutions. Of course, the solutions should not all be too similar: then there is little to gain from combining different solutions. The top-left Diversity plot tracks the average diversity of solutions in each of the feasible and infeasible solution populations. The Objectives plot gives an overview of the best and average solution quality in the current population. The bottom-left figure shows iteration runtimes in seconds. Finally, the Solution plot shows the best observed solution.

The VRP with time windows¶

Reading the instance¶

We start with a basic example that loads an instance and solves it using the standard configuration used by the solve method. For the basic example we use one of the well-known Solomon instances.

We again use the function pyvrp.read. We pass the argument round_func="dimacs" following the DIMACS VRP challenge convention, this computes distances and durations truncated to one decimal place.

[7]:

INSTANCE = read("data/RC208.vrp", round_func="dimacs")
BKS = read_solution("data/RC208.sol")

Let’s plot the instance and see what we have. The function plot_instance will plot time windows, delivery demands and coordinates, which should give us a good impression of what the instance looks like. These plots can also be produced separately by calling the appropriate plot_* function: see the API documentation for details.

[8]:

fig = plt.figure(figsize=(12, 6))
plot_instance(INSTANCE, fig)

Solving the instance¶

[9]:

result = solve(INSTANCE, stop=MaxIterations(1000), seed=42, display=False)
print(result)

Solution results
================
    # routes: 4
     # trips: 4
   # clients: 100
   objective: 7765
    distance: 7765
    duration: 17765
# iterations: 1000
    run-time: 3.19 seconds

Routes
------
Route #1: 90 82 99 52 57 49 19 18 48 21 23 25 77 58 75 97 59 87 86 74 24 22 20 66
Route #2: 65 83 64 51 85 63 76 89 33 30 32 28 26 27 29 31 34 50 62 84 56 91
Route #3: 92 95 67 71 72 41 39 38 37 35 36 40 43 44 42 61 81 54 96 93 94 80
Route #4: 69 98 88 2 6 7 79 73 78 12 14 47 17 16 15 13 9 11 10 53 60 8 46 4 45 5 3 1 70 100 55 68

[10]:

cost = result.cost() / 10
gap = 100 * (cost - BKS["cost"]) / BKS["cost"]
print(f"Found a solution with cost: {cost}.")
print(f"This is {gap:.1f}% worse than the optimal solution,", end=" ")
print(f"which is {BKS['cost']}.")

Found a solution with cost: 776.5.
This is 0.1% worse than the optimal solution, which is 776.1.

We’ve managed to find a (near) optimal solution in a few seconds!

[11]:

fig = plt.figure(figsize=(15, 9))
plot_result(result, INSTANCE, fig)
fig.tight_layout()

We can also inspect some statistics of the different routes, such as route distance, various durations, the number of stops and total delivery amount.

[12]:

solution = result.best
routes = solution.routes()

data = [
    {
        "num_stops": len(route),
        "distance": route.distance(),
        "service_duration": route.service_duration(),
        "wait_duration": route.wait_duration(),
        "time_warp": route.time_warp(),
        "delivery": route.delivery(),
    }
    for route in routes
]

tabulate(data, headers="keys", tablefmt="html")

[12]:

We can inspect the routes in more detail using the plot_route_schedule function. This will plot distance on the x-axis, and time on the y-axis, separating actual travel/driving time from waiting and service time. The clients visited are plotted as grey vertical bars indicating their time windows. In some cases, there is slack in the route indicated by a semi-transparent region on top of the earliest time line. The grey background indicates the remaining load of the truck during the route, where the (right) y-axis ends at the vehicle capacity.

[13]:

fig, axarr = plt.subplots(2, 2, figsize=(15, 9))
for idx, (ax, route) in enumerate(zip(axarr.flatten(), routes)):
    plot_route_schedule(
        INSTANCE,
        route,
        title=f"Route {idx}",
        ax=ax,
        legend=idx == 0,
    )

fig.tight_layout()

Each route begins at a given start_time, that can be obtained as follows. Note that this start time is typically not zero, that is, routes do not have to start immediately at the beginning of the time horizon.

[14]:

solution = result.best
shortest_route = min(solution.routes(), key=len)

shortest_route.start_time()

[14]:

2176

Some of the statistics presented in the plots above can also be obtained from the route schedule, as follows:

[15]:

data = [
    {
        "location": visit.location,  # Client or depot location of visit
        "start_service": visit.start_service,
        "end_service": visit.end_service,
        "service_duration": visit.service_duration,
        "wait_duration": visit.wait_duration,  # if vehicle arrives early
    }
    for visit in shortest_route.schedule()
]

tabulate(data, headers="keys", tablefmt="html")

[15]:

Solving a larger VRPTW instance¶

To show that PyVRP can also handle much larger instances, we will solve one of the largest Gehring and Homberger VRPTW benchmark instances. The selected instance - RC2_10_5 - has 1000 clients.

[16]:

INSTANCE = read("data/RC2_10_5.vrp", round_func="dimacs")
BKS = read_solution("data/RC2_10_5.sol")

[17]:

fig = plt.figure(figsize=(15, 9))
plot_instance(INSTANCE, fig)

Here, we will use a runtime-based stopping criterion: we give the solver 30 seconds to compute.

[18]:

result = solve(INSTANCE, stop=MaxRuntime(30), seed=42, display=False)

[19]:

cost = result.cost() / 10
gap = 100 * (cost - BKS["cost"]) / BKS["cost"]
print(f"Found a solution with cost: {cost}.")
print(f"This is {gap:.1f}% worse than the best-known solution,", end=" ")
print(f"which is {BKS['cost']}.")

Found a solution with cost: 27816.8.
This is 7.8% worse than the best-known solution, which is 25797.5.

[20]:

plot_result(result, INSTANCE)
plt.tight_layout()

Conclusion¶

In this notebook, we used PyVRP’s solve functionality to solve a CVRP instance with 438 clients to near-optimality, as well as several VRPTW instances, including a large 1000 client instance. Moreover, we demonstrated how to use the plotting tools to visualise the instance and statistics collected during the search procedure.