---
# System prepended metadata

title: VRP Intro
tags: [Vehicle Routing Problems, OR-Tools]

---

# VRP
We consider only the problems concerning the distribution of goods between depots and final users (customers). These problems are generally known as Vehicle Routing Problems (VRPs).

The distribution of goods concerns the service, in a given time period, of a set of customers by a set of vehicles, which are located in one or more depots, are operated by a set of drivers, and perform their movements by using an appropriate road network. In particular, the solution of a VRP calls for the determination of a set of routes, each performed by a single vehicle that starts and ends at its own depot, such that all the requirements of the customers are fulfilled, all the operational constraints are satisfied, and the global transportation cost is minimized. 

$G = (V,A)$ means a complete graph, where $V=\{0,1,...,n\}$ is the vertex set, and $A$ is the arc set. Vertex $0$ is depot.
A nonnegative $c_{ij}, (i,j)\in A$, travel cost from vertex $i$ to $j$
### Variations of VRP
There are several VRP variants that adapt to different scenario requirements:

1. **CVRP (Capacitated VRP)**: Optimizes routes when vehicles have capacity constraints.
2. **VRPTW (VRP with Time Windows)**: Includes specific time windows for each customer; vehicles must arrive within the designated time frame.
3. **MDVRP (Multi-Depot VRP)**: Optimizes routes when there are multiple depots; each vehicle must start from one depot.
4. **SDVRP (Split Delivery VRP)**: Allows splitting a customer's demand across multiple vehicles.
5. **VRPPD (Pickup and Delivery VRP)**: Includes both pickup and delivery tasks, where vehicles must follow a specified sequence.


# The Capacitated Vehicle Routing Problem (CVRP)

## Statement
- $0$ is the **depot** point.
- $N = \{1, 2,..., n\}$ denoted the **customers**.
- $q_i$ is **weight** of the goods that has to be delivered to customer $i \in N$,and $q_0 := 0$
- $K = \{1, 2,...,|K|\}$ is the **fleet**
- $|K|$ vehicles are available at the depot, have the same capacity $Q > 0$, and operating at identical costs.
- The customer subset $S \subseteq N$ denoted a vehicle **services** starts from $0$, returns to $0$, moves once to each $i$ in $S$. 
- $c_{ij}$ denoted the **travel cost** from $i$ to $j$.
- $V = \{0\} \cup N = \{0,1,..., n\}$ is the set of **vertices**.

### Underlying graph
The given information can be structured using an undirected or directed graph with the size $|K|$ of the fleet $K$ and the capacity $Q$.

- The **symmetric** case (SCVRP) 
	- Symmetric costs $c_{ij} = c_{ji}$
	- Complete and undirected weighted graph $G = (V,E,c_{ij} , q_i)$
	- Edge set $E = \{e = \{i, j \} = \{ j ,i\} | i, j \in
V,i \neq j \}$
	- Edge costs $c_{ij}$ for $\{i, j \} \in E$. 
	- $|E| = \frac{n(n + 1)}{2}$
- The **asymmetric** case (ACVRP) 
	- Asymmetric costs $c_{ij} \neq c_{ji}$
	- Complete digraph $G = (V,A,c_{ij} , q_i)$
	- Arc set $A = \{(i, j) \in V \times V | i \neq j \}$
	- Arc costs $c_{ij}$ for $(i, j) \in A$. 
	- $|A| = n(n + 1)$ 
- Both graphs contain $\mathscr{O}(n^2)$ link.

### Feasible
To a CVRP consists of $|K|$ feasible routes, one for each vehicle $k \in K$
- A **route** is a sequence $r = (i_{0},i_{1},i_{2},...,i_{s},i_{s+1})$ with $i_{0} = i_{s+1} = 0$, in which the set $S = \{i_1,...,i_s\} \subseteq N$ of customers is visited.
- The **cost of route** $c(r)=\sum_{p=0}^{s} c_{i_p i_{p+1}}$
- The route called **feasible** if 
	1. $q(S) := \sum_{i\in S}q_i \leq Q$
	2. $i_j \neq i_k$ for all $1 \leq j < k \leq s$ (no customer is visited more than once)
- such subset $S \subseteq N$ is called **feasible cluster**
- **feasible solution** if the routes $r_1, r_2,..., r_{|K|}$ are feasible and corresponding clusters $S_1, S_2, ..., S_{|K|}$ form a partition of $N$ .

Consists of two interdependent tasks:
1. The partitioning of the customer set $N$ into feasible clusters $S_1 ,..., S_{|K|}$
2. The routing of each vehicle $k \in K$ through $\{0\} \cup S_k$




## Important mathematical programming formulations

### Basic Notation
Let $S \subseteq V$ be an arbitrary subset of vertices. 
- For undirected graphs, the **cut** set $\delta(S) = \{\{i, j \} \in E | i \in S, j \in S\}$ (set $E(S) = \{\{i, j\} \in E | i, j \in S\}$) is the set of edges with exactly one (both) endpoint(s) in $S$. 
- For directed graphs $G = (V,A)$, 
	- **in-arcs** of $S$ is $\delta^−(S) = \{(i, j) \in A | i \not \in S, j ∈ S\}$
	- **out-arcs** is $\delta^+(S) = \{(i, j) \in A |i \in S, j \not \in S\}$. 
- $\delta(i) := \delta(\{i\})$ for singleton sets $S = \{i\}$ (similarly, $\delta^+(i)$ and $\delta^−(i)$). 
- $A(S) = \{(i, j) \in A | i, j \in S\}$ is the set of all arcs connecting vertices in $S$
- For any vector of variables or coefficients $x$ indexed by $i \in J$ and any subset $I \subseteq J$, the term $x(I)$ means $x(I) = \sum_{i\in I} x_i$

For a customer subset $S \subseteq N$, 
- let $r(S)$ be the minimum number of vehicle routes needed to serve $S$. which can be computed by solving a **binpacking problem** with items $N$ of weight $q_i$ ,$i \in N$, and bins of size $Q$. 
- A **lower bound** is given by $⌈\frac{q(S)}{Q}⌉$ .

### Compact Formulations
#### The two-index vehicle-flow formulations:
- $x_{ij}$ is the movement of the vehicle between i and j (from i to j) for $\{i, j\} \in E$ (or $(i, j) \in A$)
1. For the directed CVRP
	:::	success
	 **VRP1**
	The minimization of the overall routing costs:
	 $$\min\sum_{(i,j) \in A} c_{ij}x_{ij}$$
	 s.t.
	 \begin{align}
	 \begin{aligned}
		x(\delta^+(i))&:=\sum_{i \in \delta^+(i)}x_{ij}=1 && \forall i \in N, \\
		x(\delta^-(j))&:=\sum_{j \in \delta^-(j)}x_{ij}=1 &&  \forall j \in N, \\
		x(\delta^+(0))&:=\sum_{j \in \delta^+(0)}x_{0j}=|K|, && \\
	   x(\delta^+(S)) &:=\sum_{(i,j) \in \delta^+(S)}x_{ij}\ge r(S) && \forall S \subseteq N, S\not = \emptyset, \\
		& x_{ij}\in \{0,1\}  && \forall (i,j) \in A, 
	 \end{aligned}
	 \end{align}

	:::
2. For the undirected CVRP
	:::	success
	**VRP2**
	The minimization of the overall routing costs:
	 $$\min\sum_{(i,j) \in E} c_{ij}x_{ij}$$
	 s.t.
	 \begin{align}
	 \begin{aligned}
		x(\delta(i))&:=\sum_{i \in \delta(i)}x_{ij}=2 && \forall i \in N, \\
		x(\delta(0))&:=\sum_{j \in \delta(0)}x_{0j}=2|K|, && \\
	   x(\delta(S)) &:=\sum_{(i,j) \in \delta(S)}x_{ij}\ge 2r(S) && \forall S \subseteq N, S\not = \emptyset, \\
		& x_{ij}\in \{0,1,2\}  && \forall (i,j) \in \delta(0), \\
		& x_{ij}\in \{0,1\}  && \forall (i,j) \in E \backslash \delta(0), 
	 \end{aligned}
	 \end{align}

	:::

#### The three-index vehicle-flow formulation:
- based on a directed graph G = (V,A)
- replace the depot point $0$ by $o$ (**startpoint**) and $d$ (**endpoint**) of a route
- the vertex sets is $V := N \cup \{o, d\}$
- the arc sets $A := (V \backslash \{d\}) × (V \backslash \{o\})$
- $x_{ijk}$ show the movement of the vehicles over the arcs
- $x_{ijk}=1$ iff vehicle $k$ moves over the arc $(i, j) \in A$
- $y_{ik}=1$ iff vehicle $k$ visits the vertex $i \in V$
- $u_{ik}$ shows the load in vehicle $k$ directly before reaching $i \in N$
- $q_0 = q_d = 0$
:::	success
**VRP3**
The minimization of the overall routing costs:
$$\min\sum_{k \in K}\sum_{(i,j) \in A} c_{ij}x_{ijk}$$
s.t.
\begin{align}
\begin{aligned}
& \sum_{k \in K}y_{ik}=1 && \forall i \in N\\
x_k(\delta^+(i))-x_k(\delta^-(i))&=\sum_{i \in \delta^+(i)}x_{ijk}-\sum_{i \in \delta^-(i)}x_{ijk}, &&\\
&=\begin{cases}
1, &i=o\\
0, &i\in N
\end{cases} && \forall i \in V\backslash\{d\},k \in K, \\
y_{ik}&=x_k(\delta^+(i))&& \forall i \in V\backslash\{d\},k \in K, \\
y_{dk}&=x_k(\delta^-(d))&& \forall k \in K, \\
u_{ik}-u_{jk}+Qx_{ijk}&\leq Q-q_j&& \forall (i,j) \in A,k \in K, \\
q_i&\leq u_{ik} \leq Q&&\forall i \in V,k \in K, \\
x&:= (x_{k})\in \{0,1\}^{K \times A}, \\
y&:= (y_{k})\in \{0,1\}^{K \times V}  
\end{aligned}
\end{align}
:::



# Vehicle Routing Problem with Time Windows, VRPTW
## Introduction
The Vehicle Routing Problem with Time Windows (VRPTW) is an extension of the Capacitated Vehicle Routing Problem (CVRP) where each customer must be served within a specific time interval, called a time window. This problem arises in various industries like security patrol, bank deliveries, postal deliveries, refuse collection, grocery delivery, school bus routing, and newspaper distribution. If all time windows are very wide, the VRPTW becomes similar to the CVRP.
### Challenges from VRPTW
- **Much more complicated than CVRP**
  In CVRP, it is easy to find the smallest routes number you have(${\lceil}\sum_{i\in S} q_i/Q{\rceil}$), but in VRPTW this is non-trivial.
- 
## Mathematical Formulations

\begin{aligned}
\min &\sum_{k\in K}\sum_{(i,j)\in A} c_{ij}x_{ijk} &&  \\
\text{s.t.} &\sum_{k \in K}\sum_{j\in \delta^+(i)} x_{ijk} = 1, && \forall i\in N, k\in K && \\
&\sum_{j\in \delta^+(0)} x_{0jk} = 1, && \forall k\in K && \\
&\sum_{i\in \delta^-(j)} x_{ijk} - \sum_{i\in \delta^+(j)} x_{jik} = 0, && \forall k\in K, j\in N && \\
&\sum_{i\in \delta^-(n+1)} x_{i,n+1,k} = 1, && \forall k\in K && \\
&x_{ijk}(T_{ik} + s_i + t_{ij} - T_{jk}) \leq 0, && \forall k\in K, (i,j)\in A && \\
&a_i \leq T_{ik} \leq b_i, && \forall k\in K, i\in V && \\
&\sum_{i \in N}q_i\sum_{j\in \delta^+(i)} x_{ijk} \leq Q, && \forall k\in K, i\in N && \\
&x_{ijk} \in \{0,1\}, && \forall k\in K, (i,j)\in A &&
\end{aligned}
where:
- $T_{ik}$ is time variable representing the start of service at vertex i by vehicle k.
- $a_i$ and $b_i$ are the early and late start times for service at vertex i, respectively.
- $s_{i}$ is the service time at vertex i.
- $t_{ij}$ is the travel time from vertex i to vertex j.

# Solving VRP with OR-Tools


## Download the necessary libraries


```python=
!pip install osmnx
!pip install ortools
```


## Import libraries


```python=
import networkx as nx
import osmnx as ox
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from copy import deepcopy
from typing import Tuple, List, Dict, Any, Callable
```

## Set global parameters and data preprocessing
### Map



```python=+
# Download road network
map = ox.graph_from_place("竹中里, 前鎮區, 高雄市, 台灣", network_type="all")

# Plot the original map
fig, ax = plt.subplots(figsize=(10, 10))
ox.plot_graph(map, ax=ax, node_color='black', node_size=10, bgcolor='white', show=False, close=False)

# Select some nodes as nodes in the VRP problem and create a subgraph
nodes = list(map.nodes)[:10]  # Select the first 10 nodes for testing

# Get the xy coordinates of the selected nodes
node_xs = [map.nodes[node]['x'] for node in nodes]
node_ys = [map.nodes[node]['y'] for node in nodes]

# Mark the selected nodes in red
ax.scatter(node_xs, node_ys, c='red', s=50, zorder=5)

# Show the map
plt.show()
```


    
![Map](https://i.imgur.com/mvWbNbV.png)
    


Select the first ten points of the node sequence in the graph as customers



```python=+
# Default road types and their speeds (unit: km/h)
hwy_speeds = {
    'residential': 30,
    'primary': 60,
    'secondary': 50,
    'tertiary': 40,
    'motorway': 100
}

# Set speed limit, compute expected travel time
map = ox.add_edge_speeds(map, hwy_speeds=hwy_speeds)
map = ox.add_edge_travel_times(map)

# Convert into a DiGraph
map_di = nx.DiGraph(map)
for u, v, k, data in map.edges(keys=True, data=True):
    if not map_di.has_edge(u, v) or data['travel_time'] < map_di[u][v]['travel_time']:
        map_di.add_edge(u, v, **data)
```

### Distance Matrix


```python=+
# Calculate the shortest path distance matrix between nodes
distance_matrix = nx.floyd_warshall_numpy(map, nodelist=map.nodes)
# Intercept the distance matrix of the first ten points of the node sequence
distance_matrix = distance_matrix[:10,:10]
```

## Use Ortools to solve several types of VRP

### Custom Libraries

#### OR-Tools related


```python=+
from ortools.constraint_solver import pywrapcp, routing_enums_pb2

def create_data_model(distance_matrix, vrp_type = None):
    """Stores the data for the problem."""
    match vrp_type:
        case "vrp":
            data = {}
            data['distance_matrix'] = distance_matrix.tolist()  # Convert numpy matrix to list format
            data['num_vehicles'] = 4  # Assuming there are 4 vehicles
            data['depot'] = 0  # Assuming the starting point is the first node
            return data

        case "cvrp":
            data = {}
            data["distance_matrix"] = distance_matrix.tolist()  # Convert numpy matrix to list format
            data["demands"] = [0, 6, 3, 2, 1, 2, 4, 4, 8, 8]
            data["vehicle_capacities"] = [9, 9, 10, 10]
            data["num_vehicles"] = 4
            data["depot"] = 0
            return data

        case "vrptw":
            data = {}
            data["time_matrix"] = distance_matrix.tolist()  # Convert numpy matrix to list format
            data["time_windows"] = [
                (0, 5),  # depot
                (7, 12),  # 1
                (10, 15),  # 2
                (16, 18),  # 3
                (10, 13),  # 4
                (0, 5),  # 5
                (5, 10),  # 6
                (0, 4),  # 7
                (5, 10),  # 8
                (0, 3),  # 9
            ]
            data["num_vehicles"] = 4
            data["depot"] = 0
            return data

        case _:
            print("You should input the type of the problem: {\"vrp\", \"cvrp\", \"vrptw\"}")

def print_solution(manager, routing, solution):
    """Return solution on console and returns the route."""
    routes = []
    for vehicle_id in range(manager.GetNumberOfVehicles()):
        index = routing.Start(vehicle_id)
        route = []
        while not routing.IsEnd(index):
            route.append(manager.IndexToNode(index))
            index = solution.Value(routing.NextVar(index))
        route.append(manager.IndexToNode(index))
        routes.append(route)
    return routes
```

#### Define the function for map visualization


```python=+
def deterministic_shortest_path(graph: nx.DiGraph, origin: int, destination: int) -> Tuple[List[int], float]:
    """
    Compute the deterministic shortest path between origin and destination.

    Args:
    graph (nx.DiGraph): The road network graph.
    origin (int): The starting node.
    destination (int): The ending node.

    Returns:
    Tuple[List[int], float]: The shortest path and its total travel time.
    """
    if origin not in graph or destination not in graph:
        raise ValueError("Origin or destination not in graph.")

    # Compute value function without uncertainty using Bellman-Ford algorithm
    value_dict = {node: float('inf') for node in graph.nodes()}
    value_dict[destination] = 0

    for _ in range(len(graph) - 1):
        for u, v, data in graph.edges(data=True):
            if value_dict[u] > value_dict[v] + data['travel_time']:
                value_dict[u] = value_dict[v] + data['travel_time']

    # Compute shortest path by value function
    shortest_path = [origin]
    shortest_path_length = 0.0

    while shortest_path[-1] != destination:
        current = shortest_path[-1]
        next_node = min(
            graph.successors(current),
            key=lambda x: graph[current][x]['travel_time'] + value_dict[x]
        )
        shortest_path.append(next_node)
        shortest_path_length += graph[current][next_node]['travel_time']

    return shortest_path, shortest_path_length


def visualize_path(graph: nx.DiGraph, path: List[int]) -> None:
    """
    Visualize the shortest path on the graph.

    Args:
    graph (nx.DiGraph): The road network graph.
    path (List[int]): The shortest path to visualize.
    """
    fig, ax = ox.plot_graph_route(
        nx.MultiDiGraph(graph),
        path,
        route_linewidth=2,
        node_size=7,
        bgcolor='w',
        node_color='k',
        edge_color='#696969'
    )
    plt.show()

def visualize_paths(graph: nx.DiGraph, paths: List[int]) -> None:
    """
    Visualize the shortest path on the graph.

    Args:
    graph (nx.DiGraph): The road network graph.
    path (List[int]): The shortest path to visualize.
    """
    fig, ax = ox.plot_graph_routes(
        nx.MultiDiGraph(graph),
        paths,
        route_linewidth=2,
        node_size=7,
        bgcolor='w',
        node_color='k',
        edge_color='#696969'
    )
    plt.show()
	
def deterministic_shortest_route(graph: nx.DiGraph, route: List[int]) -> Tuple[List[List[int]], List[float]]:
    '''
    Calculate the shortest path and time cost of each section
    '''
    shortest_paths = []
    travel_times = []
    for i in range(len(route)-1):
        shortest_path, travel_time = deterministic_shortest_path(graph, route[i], route[i+1])
        shortest_paths.append(shortest_path)
        travel_times.append(travel_time)

    return shortest_paths, travel_times
```

### Simple VRP


```python=+
def solve_vrp(distance_matrix):
    """Solves the VRP problem."""
    # Create the data model
    data = create_data_model(distance_matrix, vrp_type="vrp")

    # Create the routing index manager
    manager = pywrapcp.RoutingIndexManager(len(data['distance_matrix']),
                                           data['num_vehicles'], data['depot'])

    # Create Routing Model
    routing = pywrapcp.RoutingModel(manager)

    # Create and register a transit callback
    def distance_callback(from_index, to_index):
        """Returns the distance between the two nodes."""
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return data['distance_matrix'][from_node][to_node]

    transit_callback_index = routing.RegisterTransitCallback(distance_callback)

    # Define cost of each arc
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

    # Add Distance constraint
    dimension_name = 'Distance'
    routing.AddDimension(
        transit_callback_index,
        0,  # no slack
        3000,  # maximum distance per vehicle
        True,  # start cumul to zero
        dimension_name)
    distance_dimension = routing.GetDimensionOrDie(dimension_name)
    distance_dimension.SetGlobalSpanCostCoefficient(100)

    # Set first solution heuristic
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)

    # Solve the problem
    solution = routing.SolveWithParameters(search_parameters)

    # return solution
    if solution:
        routes = print_solution(manager, routing, solution)
        return routes
    else:
        print("No solution found!")
        return None
```


```python=+
# Solve VRP problem
routes = solve_vrp(distance_matrix)
# Replace the values in routes with the values on the map
routes = [[nodes[i] for i in route] for route in routes]
```


```python=+
def graph_ordered_nodes(map, selected_nodes: List[int], paths: List[int]):
    # Get the location of nodes and edges
    nodes, edges = ox.graph_to_gdfs(map)

    # Create a map
    fig, ax = ox.plot_graph_routes(
                            map,
                            paths,
                            show=False,
                            close=False,
                            route_linewidth=2,
                            node_size=7,
                            bgcolor='w',
                            node_color='k',
                            edge_color='#696969')

    # Add number labels to the selected nodes, starting from 0, excluding the last one
    for i, node_id in enumerate(selected_nodes[:-1]):
        if node_id in map.nodes:
            data = map.nodes[node_id]
            x, y = data['x'], data['y']
            ax.text(x, y, str(i), fontsize=12, ha='center', va='center',
                    bbox=dict(facecolor='white', edgecolor='none', alpha=0.7))
        else:
            print(f"Node {node_id} not found in the graph.")

    # Adjust layout to fit legend
    plt.tight_layout()
    plt.subplots_adjust(right=0.8)

    plt.show()
```

#### Visualization


```python=+
shortest_paths, travel_times = deterministic_shortest_route(map_di,routes[3])
graph_ordered_nodes(map, routes[3], shortest_paths)
```


    
![VRP](https://i.imgur.com/6fO2GcK.png)
    


### CVRP


```python=+
def solve_cvrp(distance_matrix):
    """Solve the CVRP problem."""
    # Instantiate the data problem.
    data = create_data_model(distance_matrix,vrp_type="cvrp")

    # Create the routing index manager.
    manager = pywrapcp.RoutingIndexManager(
        len(data["distance_matrix"]), data["num_vehicles"], data["depot"]
    )

    # Create Routing Model.
    routing = pywrapcp.RoutingModel(manager)

    # Create and register a transit callback.
    def distance_callback(from_index, to_index):
        """Returns the distance between the two nodes."""
        # Convert from routing variable Index to distance matrix NodeIndex.
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return data["distance_matrix"][from_node][to_node]

    transit_callback_index = routing.RegisterTransitCallback(distance_callback)

    # Define cost of each arc.
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

    # Add Capacity constraint.
    def demand_callback(from_index):
        """Returns the demand of the node."""
        # Convert from routing variable Index to demands NodeIndex.
        from_node = manager.IndexToNode(from_index)
        return data["demands"][from_node]

    demand_callback_index = routing.RegisterUnaryTransitCallback(demand_callback)
    routing.AddDimensionWithVehicleCapacity(
        demand_callback_index,
        0,  # null capacity slack
        data["vehicle_capacities"],  # vehicle maximum capacities
        True,  # start cumul to zero
        "Capacity",
    )

    # Setting first solution heuristic.
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC
    )
    search_parameters.local_search_metaheuristic = (
        routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH
    )
    search_parameters.time_limit.FromSeconds(1)

    # Solve the problem.
    solution = routing.SolveWithParameters(search_parameters)

    # return solution
    if solution:
        routes = print_solution(manager, routing, solution)
        return routes
    else:
        print("No solution found!")
        return None
```


```python=+
# Solve CVRP problem
routes = solve_cvrp(distance_matrix)
routes = [[nodes[i] for i in route] for route in routes]
```

#### Visualization


```python=+
i = 0
shortest_paths, travel_times = deterministic_shortest_route(map_di,routes[i])
graph_ordered_nodes(map, routes[i], shortest_paths)
```


    
![CVRP_0](https://i.imgur.com/ebW6Sxu.png)
    



```python=+
i = 1
shortest_paths, travel_times = deterministic_shortest_route(map_di,routes[i])
graph_ordered_nodes(map, routes[i], shortest_paths)
```


    

![CVRP_1](https://i.imgur.com/mPjaYho.png)



```python=+
i = 2
shortest_paths, travel_times = deterministic_shortest_route(map_di,routes[i])
graph_ordered_nodes(map, routes[i], shortest_paths)
```


    

![CVRP_2](https://i.imgur.com/uVYKOnb.png)
    



```python=+
i = 3
shortest_paths, travel_times = deterministic_shortest_route(map_di,routes[i])
graph_ordered_nodes(map, routes[i], shortest_paths)
```

![CVRP_3](https://i.imgur.com/33bLC6K.png)
    


### VRPTW


```python=+
def solve_vrptw(distance_matrix):
    """Solve the VRP with time windows."""
    # Instantiate the data problem.
    data = create_data_model(distance_matrix,vrp_type="vrptw")

    # Create the routing index manager.
    manager = pywrapcp.RoutingIndexManager(
        len(data["time_matrix"]), data["num_vehicles"], data["depot"]
    )

    # Create Routing Model.
    routing = pywrapcp.RoutingModel(manager)

    # Create and register a transit callback.
    def time_callback(from_index, to_index):
        """Returns the travel time between the two nodes."""
        # Convert from routing variable Index to time matrix NodeIndex.
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return data["time_matrix"][from_node][to_node]

    transit_callback_index = routing.RegisterTransitCallback(time_callback)

    # Define cost of each arc.
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

    # Add Time Windows constraint.
    time = "Time"
    routing.AddDimension(
        transit_callback_index,
        30,  # allow waiting time
        30,  # maximum time per vehicle
        False,  # Don't force start cumul to zero.
        time,
    )
    time_dimension = routing.GetDimensionOrDie(time)
    # Add time window constraints for each location except depot.
    for location_idx, time_window in enumerate(data["time_windows"]):
        if location_idx == data["depot"]:
            continue
        index = manager.NodeToIndex(location_idx)
        time_dimension.CumulVar(index).SetRange(time_window[0], time_window[1])
    # Add time window constraints for each vehicle start node.
    depot_idx = data["depot"]
    for vehicle_id in range(data["num_vehicles"]):
        index = routing.Start(vehicle_id)
        time_dimension.CumulVar(index).SetRange(
            data["time_windows"][depot_idx][0], data["time_windows"][depot_idx][1]
        )

    # Instantiate route start and end times to produce feasible times.
    for i in range(data["num_vehicles"]):
        routing.AddVariableMinimizedByFinalizer(
            time_dimension.CumulVar(routing.Start(i))
        )
        routing.AddVariableMinimizedByFinalizer(time_dimension.CumulVar(routing.End(i)))

    # Setting first solution heuristic.
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC
    )

    # Solve the problem.
    solution = routing.SolveWithParameters(search_parameters)

    # return solution
    if solution:
        routes = print_solution(manager, routing, solution)
        return routes
    else:
        print("No solution found!")
        return None
```


```python=+
# Solve VRP problem
routes = solve_vrptw(distance_matrix)
routes = [[nodes[i] for i in route] for route in routes]
```

#### Visualization


```python=+
i = 2
shortest_paths, travel_times = deterministic_shortest_route(map_di,routes[i])
graph_ordered_nodes(map, routes[i], shortest_paths)
```


    
![VRPTW_2](https://i.imgur.com/iaS82UF.png)
    



```python=+
i = 3
shortest_paths, travel_times = deterministic_shortest_route(map_di,routes[i])
graph_ordered_nodes(map, routes[i], shortest_paths)
```


    
![VRPTW_3](https://i.imgur.com/HB2QPCJ.png)
