Skip to content

Capacitated Vehicle Routing Problem (CVRP) Example

Download Notebook


The Capacitated Vehicle Routing Problem asks for a set of vehicle routes, each starting and ending at a single depot, that together visit every customer exactly once without exceeding the per-vehicle capacity, minimizing total travelled distance. It is NP-hard and central to logistics and distribution planning.

import getpass
import os

import numpy as np
from dotenv import load_dotenv
from luna_quantum.algorithms import SCIP

from luna_usecases.vehicle_routing_problem import (
    CvrpCollection,
    CvrpData,
    CvrpFormulation,
    CvrpInstance,
)

load_dotenv()
if "LUNA_API_KEY" not in os.environ:
    os.environ["LUNA_API_KEY"] = getpass.getpass("Enter your Luna API key: ")

Create Data

Build an instance from 2D coordinates using the from_coordinates factory (Euclidean distances). Node 0 is the depot; demands are per customer.

coords = np.array(
    [
        [50.0, 50.0],  # depot
        [20.0, 30.0],
        [80.0, 20.0],
        [70.0, 80.0],
        [30.0, 70.0],
    ]
)
demands = [0.0, 2.0, 3.0, 2.0, 3.0]
data = CvrpData.from_coordinates(
    coords=coords,
    demands=demands,
    n_vehicles=2,
    vehicle_capacity=6.0,
)
print(data.to_string())
Capacitated Vehicle Routing Data:
  Customers: 4
  Vehicles: 2
  Vehicle capacity: 6.0
  Depot: 0
  Total demand: 10.0
  Distance matrix shape: (5, 5)

Plot Data

Visualize the depot and customer locations.

data.plot()

<Axes: title={'center': 'CVRP — 4 customers, 2 vehicles'}>
png

Create Formulation

Minimize total distance using the Miller-Tucker-Zemlin (MTZ) formulation with capacity and subtour-elimination constraints.

formulation = CvrpFormulation()
print(formulation.to_string(data))
Capacitated Vehicle Routing Problem Formulation (MTZ):
  Customers: 4
  Vehicles: 2
  Vehicle capacity: 6.0

Decision Variables:
  x[i,j] in {0,1} for i != j, i,j = 0, ..., 4
  x[i,j] = 1 if a vehicle travels from node i to node j
  u[i] integer in [demand_i, 6.0] for i = 1, ..., 4
  u[i] = cumulative load just after visiting customer i

Objective:
  minimize sum_{i!=j} d[i,j] * x[i,j]

Constraints:
  1. Each customer entered once (4 constraints):
     sum_{j!=i} x[j,i] == 1  for all customers i
  2. Each customer left once (4 constraints):
     sum_{j!=i} x[i,j] == 1  for all customers i
  3. Depot out-degree:
     sum_{j>=1} x[0,j] == 2
  4. Depot in-degree:
     sum_{j>=1} x[j,0] == 2
  5. MTZ capacity / subtour elimination:
     u[i] - u[j] + 6.0 * x[i,j] <= 6.0 - demand_j

Create Instance

Combine data and formulation into a solvable instance.

instance = CvrpInstance(data=data, formulation=formulation)
print(instance.to_string())
Data:Capacitated Vehicle Routing Data:
  Customers: 4
  Vehicles: 2
  Vehicle capacity: 6.0
  Depot: 0
  Total demand: 10.0
  Distance matrix shape: (5, 5)
Formulation:Capacitated Vehicle Routing Problem Formulation (MTZ):
  Customers: 4
  Vehicles: 2
  Vehicle capacity: 6.0

Decision Variables:
  x[i,j] in {0,1} for i != j, i,j = 0, ..., 4
  x[i,j] = 1 if a vehicle travels from node i to node j
  u[i] integer in [demand_i, 6.0] for i = 1, ..., 4
  u[i] = cumulative load just after visiting customer i

Objective:
  minimize sum_{i!=j} d[i,j] * x[i,j]

Constraints:
  1. Each customer entered once (4 constraints):
     sum_{j!=i} x[j,i] == 1  for all customers i
  2. Each customer left once (4 constraints):
     sum_{j!=i} x[i,j] == 1  for all customers i
  3. Depot out-degree:
     sum_{j>=1} x[0,j] == 2
  4. Depot in-degree:
     sum_{j>=1} x[j,0] == 2
  5. MTZ capacity / subtour elimination:
     u[i] - u[j] + 6.0 * x[i,j] <= 6.0 - demand_j

Formulate Model

Translate the instance into a mixed-integer optimization model.

model = instance.formulate()

Solve and Interpret

Solve the model with SCIP and interpret the raw result into a use-case-specific solution.

scip = SCIP()
job = scip.run(model)
sol = job.result()
uc_solution = instance.interpret(sol)
print(uc_solution.to_string())
Waiting for result | elapsed 0:00:00 | poll 1 | next check in 0.5s
Waiting for result | elapsed 0:00:01 | poll 2 | next check in 1s
Waiting for result | elapsed 0:00:03 | poll 3 | next check in 2s
Waiting for result | elapsed 0:00:05 | poll 4 | next check in 4s
Waiting for result | elapsed 0:00:09 | poll 5 | next check in 8s
Waiting for result | elapsed 0:00:18 | poll 6 | next check in 16s
2026-08-20 16:18:57 INFO     Solve job '7f0b8169-4b7d-4e7b-a4f9-b822677cde61' finished with status 'DONE' after    
                             0:00:35.                                                                              
Capacitated Vehicle Routing Solution:
  Routes (2):
    Vehicle 0: 0 -> 1 -> 2 -> 0
    Vehicle 1: 0 -> 3 -> 4 -> 0
  Total distance: 244.88
  Valid: True

Plot Solution

Visualize the vehicle routes.

uc_solution.plot(data)

<Axes: title={'center': 'CVRP Solution — distance: 244.9, valid=True'}>
png

Collections

Generate a benchmark collection of random instances for batch processing.

collection = CvrpCollection.from_random(min_customers=3, max_customers=4, num_instances=2, seed=42)
model = collection.instances[0].formulate()
print(model)
Model: vehicle_routing_problem<s42_c3_i0>
Minimize
  44.92191775374967 * x_0_1 + 44.44282935787485 * x_0_2
  + 65.2734432668789 * x_0_3 + 44.92191775374967 * x_1_0
  + 84.13995813993871 * x_1_2 + 20.71535492454286 * x_1_3
  + 44.44282935787485 * x_2_0 + 84.13995813993871 * x_2_1
  + 104.74544502819062 * x_2_3 + 65.2734432668789 * x_3_0
  + 20.71535492454286 * x_3_1 + 104.74544502819062 * x_3_2
Subject To
  in_degree_1: x_0_1 + x_2_1 + x_3_1 == 1
  in_degree_2: x_0_2 + x_1_2 + x_3_2 == 1
  in_degree_3: x_0_3 + x_1_3 + x_2_3 == 1
  out_degree_1: x_1_0 + x_1_2 + x_1_3 == 1
  out_degree_2: x_2_0 + x_2_1 + x_2_3 == 1
  out_degree_3: x_3_0 + x_3_1 + x_3_2 == 1
  depot_out_degree: x_0_1 + x_0_2 + x_0_3 == 2
  depot_in_degree: x_1_0 + x_2_0 + x_3_0 == 2
  mtz_1_2: 9 * x_1_2 + u_1 - u_2 <= 5
  mtz_1_3: 9 * x_1_3 + u_1 - u_3 <= 6
  mtz_2_1: 9 * x_2_1 - u_1 + u_2 <= 6
  mtz_2_3: 9 * x_2_3 + u_2 - u_3 <= 6
  mtz_3_1: 9 * x_3_1 - u_1 + u_3 <= 6
  mtz_3_2: 9 * x_3_2 - u_2 + u_3 <= 5
Bounds
  3 <= u_1 <= 9
  4 <= u_2 <= 9
  3 <= u_3 <= 9
Binary
  x_0_1 x_0_2 x_0_3 x_1_0 x_1_2 x_1_3 x_2_0 x_2_1 x_2_3 x_3_0 x_3_1 x_3_2
Integer
  u_1 u_2 u_3