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 Model, read
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()
../_images/examples_basic_vrps_6_0.png

Solving the instance

We will again use the Model interface to solve the instance. The Model interface supports a convenient from_data method that can be used to instantiate a model from a known ProblemData object.

[4]:

model = Model.from_data(INSTANCE)
result = model.solve(stop=MaxIterations(2000), seed=42)
print(result)

Solution results
================
    # routes: 37
   # clients: 438
   objective: 36802.00
# iterations: 2000
    run-time: 583.42 seconds

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


[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: 36802.
This is 1.1% 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()
../_images/examples_basic_vrps_11_0.png

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 Model interface. For the basic example we use one of the well-known Solomon instances.

We again use the function pyvrp.read. We pass the following arguments:

  • instance_format="solomon": this parses the instance file as a Solomon formatted instance.

  • round_func="trunc1": following the DIMACS VRP challenge convention, this computes distances and durations truncated to one decimal place.

[7]:

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

Let’s plot the instance and see what we have. The function plot_instance will plot time windows, 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)
../_images/examples_basic_vrps_16_0.png

Solving the instance

We will again use the Model interface to solve the instance.

[9]:

model = Model.from_data(INSTANCE)
result = model.solve(stop=MaxIterations(1000), seed=42)
print(result)

Solution results
================
    # routes: 4
   # clients: 100
   objective: 7761.00
# iterations: 1000
    run-time: 27.48 seconds

Routes
------
Route #1: 90 65 82 99 52 83 64 49 19 18 48 21 23 25 77 58 75 97 59 87 74 86 57 24 22 20 66
Route #2: 94 92 95 67 62 50 34 31 29 27 26 28 30 32 33 76 89 63 85 51 84 56 91 80
Route #3: 61 42 44 39 38 36 35 37 40 43 41 72 71 93 96 54 81
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.1.
This is 0.0% worse than the optimal solution, which is 776.1.

We’ve managed to find the optimal solution in a few seconds!

[11]:

fig = plt.figure(figsize=(15, 9))
plot_result(result, INSTANCE, fig)
fig.tight_layout()
../_images/examples_basic_vrps_21_0.png

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

[12]:

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

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

header = list(data[0].keys())
rows = [datum.values() for datum in data]
tabulate(rows, header, tablefmt="html")

[12]:
num_stops distance service_duration wait_duration time_warp demand
27 2187 2700 0 0 465
24 1983 2400 0 0 381
17 1325 1700 0 0 286
32 2266 3200 0 0 592

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. We can see a jump to the start of the time window in the main (earliest) time line when a vehicle arrives early at a customer and has to wait. 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.reshape(-1), routes)):
    plot_route_schedule(
        INSTANCE,
        route,
        title=f"Route {idx}",
        ax=ax,
        legend=idx == 0,
    )

fig.tight_layout()
../_images/examples_basic_vrps_25_0.png

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.

[14]:

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

[15]:

fig = plt.figure(figsize=(15, 9))
plot_instance(INSTANCE, fig)
../_images/examples_basic_vrps_28_0.png

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

[16]:

model = Model.from_data(INSTANCE)
result = model.solve(stop=MaxRuntime(30), seed=42)

[17]:

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: 30963.1.
This is 20.0% worse than the best-known solution, which is 25797.5.

[18]:

plot_result(result, INSTANCE)
../_images/examples_basic_vrps_32_0.png

Conclusion

In this notebook, we used PyVRP’s Model interface 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.