Vehicle Routing Problem (VRP) Example

The Capacitated Vehicle Routing Problem (CVRP) is a combinatorial optimization problem where a fleet of vehicles with limited capacity must deliver goods to a set of customers, minimizing total travel distance.

Mathematical Formulation

Sets and indices:

  • \(V = \{0, 1, \dots, n\}\) — set of nodes, where \(0\) is the depot and \(\{1, \dots, n\}\) are customers

  • \(K = \{1, \dots, m\}\) — set of vehicles

  • \(A = \{(i, j) : i, j \in V, \; i \neq j\}\) — set of arcs

Parameters:

  • \(c_{ij}\) — travel cost (distance) from node \(i\) to node \(j\)

  • \(q_i\) — demand of customer \(i\), with \(q_0 = 0\)

  • \(Q\) — vehicle capacity

Decision variables:

\[\begin{split}x_{ijk} = \begin{cases} 1, & \text{if vehicle } k \text{ travels from node } i \text{ to node } j \\ 0, & \text{otherwise} \end{cases}\end{split}\]

Objective — minimize total travel distance:

\[\min \sum_{k \in K} \sum_{(i,j) \in A} c_{ij} \, x_{ijk}\]

Subject to:

Each customer is visited exactly once:

\[\sum_{k \in K} \sum_{j \in V} x_{ijk} = 1, \quad \forall \, i \in V \setminus \{0\}\]

Flow conservation — each vehicle that enters a node must leave it:

\[\sum_{i \in V} x_{ijk} = \sum_{i \in V} x_{jik}, \quad \forall \, j \in V, \; \forall \, k \in K\]

Each vehicle leaves and returns to the depot at most once:

\[\sum_{j \in V \setminus \{0\}} x_{0jk} \leq 1, \quad \forall \, k \in K\]

Vehicle capacity constraint:

\[\sum_{i \in V \setminus \{0\}} q_i \sum_{j \in V} x_{ijk} \leq Q, \quad \forall \, k \in K\]

Subtour elimination (for each subset \(S \subseteq V \setminus \{0\}\), \(S \neq \emptyset\)):

\[\sum_{i \in S} \sum_{j \notin S} x_{ijk} \geq 1, \quad \forall \, S \subseteq V \setminus \{0\}, \; \forall \, k \in K\]

Solution Approach

Since the CVRP is NP-hard, this implementation uses a heuristic approach:

  1. Clarke-Wright Savings Algorithm — constructive heuristic that builds an initial solution by iteratively merging routes based on the savings \(s_{ij} = c_{0i} + c_{0j} - c_{ij}\)

  2. Local search improvement — applies 2-opt (segment reversal), or-opt (chain relocation), inter-route customer relocation, and inter-route customer swap moves until no further improvement is found

[1]:
import sys
sys.path.insert(0, "../src")
[2]:
import numpy as np
from optymus.routing import solve_vrp, VRPSolver
/Users/holisticai008/Documents/repos/quantsci/optymus/.venv/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

1. Define the problem

We create an instance with 25 customers and a depot (node 0). Each customer has a demand, and each vehicle has a capacity of 30 units.

[3]:
# Depot at index 0, followed by 25 customers
np.random.seed(49)
n_customers = 25
coordinates = np.vstack([
    [50, 50],  # depot
    np.column_stack([
        np.random.randint(5, 95, size=n_customers),
        np.random.randint(5, 95, size=n_customers),
    ])
]).astype(float)

# Customer demands (depot demand = 0)
demands = np.concatenate([[0], np.random.randint(1, 10, size=n_customers)]).astype(float)

vehicle_capacity = 30

print(f"Number of customers: {n_customers}")
print(f"Total demand: {demands.sum():.0f}")
print(f"Vehicle capacity: {vehicle_capacity}")
print(f"Min vehicles needed: {int(np.ceil(demands.sum() / vehicle_capacity))}")
Number of customers: 25
Total demand: 139
Vehicle capacity: 30
Min vehicles needed: 5

2. Solve with the functional API

[4]:
result = solve_vrp(
    coordinates=coordinates,
    demands=demands,
    vehicle_capacity=vehicle_capacity,
)

result
Initial solution: 5 routes, distance = 555.3417
Final solution:   5 routes, distance = 521.3878 (6.1% improvement)
[4]:
OptimizeResult: CVRP (Clarke-Wright Savings + 2-opt)
----------------------------------------------------
  Optimal Solution   5 routes
  Function Value     5.213878e+02
  Iterations         2
  Termination        local_search_converged
  Time (s)           0.0019
  Converged          True
----------------------------------------------------

3. Inspect route details

[5]:
for i, rd in enumerate(result.route_details):
    print(f"Route {i + 1}: depot -> {rd['customers']} -> depot  "
          f"(distance={rd['distance']:.2f}, load={rd['load']:.0f}/{vehicle_capacity})")
Route 1: depot -> [6, 20, 19, 12, 17] -> depot  (distance=112.88, load=28/30)
Route 2: depot -> [23, 24, 22, 4, 21] -> depot  (distance=88.03, load=29/30)
Route 3: depot -> [8, 3, 13, 10, 16] -> depot  (distance=114.14, load=25/30)
Route 4: depot -> [14, 18, 5, 7, 11, 9] -> depot  (distance=128.23, load=28/30)
Route 5: depot -> [15, 25, 2, 1] -> depot  (distance=78.10, load=29/30)

4. Visualize the routes

[6]:
try:
    import matplotlib.pyplot as plt

    colors = plt.cm.Set2(np.linspace(0, 1, len(result.routes)))

    fig, ax = plt.subplots(figsize=(8, 8))

    # Plot depot
    ax.plot(*coordinates[0], 's', color='red', markersize=12, label='Depot', zorder=5)

    # Plot each route
    for idx, route in enumerate(result.routes):
        full_route = [0] + route + [0]
        xs = coordinates[full_route, 0]
        ys = coordinates[full_route, 1]
        ax.plot(xs, ys, 'o-', color=colors[idx], linewidth=2, markersize=8,
                label=f'Route {idx + 1} (load={result.route_details[idx]["load"]:.0f})')

    # Label customers
    for i in range(1, len(coordinates)):
        ax.annotate(f'{i}', coordinates[i], textcoords='offset points',
                    xytext=(5, 5), fontsize=9)

    ax.set_title(f'CVRP Solution — Total distance: {result.total_distance:.2f}')
    ax.legend(loc='upper left')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.show()

except ImportError:
    print("Install matplotlib for route visualization: pip install matplotlib")
Install matplotlib for route visualization: pip install matplotlib

5. Class-based API

For more control, use VRPSolver directly.

[7]:
solver = VRPSolver(
    coordinates=coordinates,
    demands=demands,
    vehicle_capacity=vehicle_capacity,
    num_vehicles=5,
    verbose=True,
)

result = solver.solve(max_iter=200)
result
Initial solution: 5 routes, distance = 555.3417
Final solution:   5 routes, distance = 521.3878 (6.1% improvement)
[7]:
OptimizeResult: CVRP (Clarke-Wright Savings + 2-opt)
----------------------------------------------------
  Optimal Solution   5 routes
  Function Value     5.213878e+02
  Iterations         2
  Termination        local_search_converged
  Time (s)           0.0012
  Converged          True
----------------------------------------------------

6. Using a distance matrix directly

If you already have a distance (or cost) matrix, pass it instead of coordinates.

[8]:
from optymus.routing import compute_distance_matrix

D = compute_distance_matrix(coordinates)

result = solve_vrp(
    distance_matrix=D,
    demands=demands,
    vehicle_capacity=vehicle_capacity,
    verbose=False,
)

print(f"Total distance: {result.total_distance:.2f}")
print(f"Vehicles used: {result.num_vehicles_used}")
print(f"Converged: {result.converged}")
Total distance: 521.39
Vehicles used: 5
Converged: True