Thana B.

Posted on Oct 6

A Practical Review: Solving Vehicle Routing Problems with OR-Tools and SCIP

#ortools

Introduction

Vehicle Routing Problems (VRP) are a class of combinatorial optimization problems that are crucial in the fields of logistics and transportation. The objective is to determine the optimal set of routes for a fleet of vehicles to traverse in order to deliver goods to a given set of customers. The challenge lies in balancing multiple constraints such as vehicle capacities, delivery time windows, and minimizing the total travel time or distance.

In this blog, we will explore how to solve VRP using two powerful optimization tools: Google OR-Tools and SCIP. Google OR-Tools is an open-source software suite for optimization, providing a range of solvers for different types of problems. SCIP, on the other hand, is a framework for Constraint Integer Programming and Mixed Integer Programming, known for its flexibility and performance.

We will implement the VRP with time windows using both OR-Tools and SCIP, compare their performance, and share insights from our experience with these tools.

Problem Statement

The specific VRP we are tackling involves a fleet of vehicles that must deliver goods to a set of customers within specified time windows. The problem can be summarized as follows:

Time Matrix: A matrix representing the travel time between each pair of locations.
Time Windows: Each location (including the depot) has a time window within which the delivery must be made.
Demands: Each customer has a demand that must be satisfied by the vehicles.
Vehicle Capacities: Each vehicle has a maximum capacity that cannot be exceeded.
Number of Vehicles: The fleet consists of a fixed number of vehicles.
Depot: The starting and ending point for all vehicles.

Note that waiting time is allowed and vehicles are not necessary to depart from its depot at 0th minute.

The objective is to minimize the total travel time while satisfying all the constraints. The problem is modeled and solved using Google OR-Tools, and the solution is compared with various models implemented in SCIP.

Tools Overview

Google OR-Tools

Brief Introduction to OR-Tools

Google OR-Tools is an open-source software suite for optimization, developed by Google. It provides a range of solvers for different types of optimization problems, including linear programming, mixed-integer programming, constraint programming, and routing problems. OR-Tools is designed to be easy to use and integrate into various applications, making it a popular choice for both academic research and industrial applications.

Key Features and Advantages

Ease of Use: OR-Tools provides a high-level API that simplifies the process of defining and solving optimization problems. This makes it accessible to users with varying levels of expertise in optimization.
Versatility: It supports a wide range of problem types, including VRP, scheduling, and resource allocation.
Performance: OR-Tools is optimized for performance, leveraging advanced algorithms and heuristics to find solutions quickly.
Cross-Platform: It is available on multiple platforms, including Windows, macOS, and Linux.
Integration: OR-Tools can be easily integrated with other Google services and tools, such as Google Cloud Platform.

Use Cases and Applications

Logistics and Transportation: Solving VRP, optimizing delivery routes, and managing fleet operations.
Scheduling: Creating efficient schedules for employees, machines, or other resources.
Resource Allocation: Optimizing the allocation of resources in various industries, such as manufacturing and healthcare.
Research and Education: Used in academic research and teaching to demonstrate optimization techniques and solve complex problems.

SCIP

Brief Introduction to SCIP

SCIP (Solving Constraint Integer Programs) is a framework for Constraint Integer Programming and Mixed Integer Programming. Developed by the Zuse Institute Berlin, SCIP is known for its flexibility and high performance. It combines techniques from constraint programming, mixed-integer programming, and SAT solving, making it a powerful tool for solving a wide range of optimization problems.

Key Features and Advantages

Flexibility: SCIP allows users to define custom constraints and objective functions, providing a high degree of control over the optimization process.
Performance: SCIP is one of the fastest non-commercial solvers available, often outperforming commercial solvers in benchmark tests.
Extensibility: Users can extend SCIP by adding custom plugins for specific problem types or optimization techniques.
Open Source: SCIP is open-source and freely available for academic and non-commercial use, making it accessible to a wide audience.
Comprehensive Documentation: SCIP provides extensive documentation and examples, helping users get started quickly and effectively.

Use Cases and Applications

Combinatorial Optimization: Solving complex combinatorial problems, such as VRP, scheduling, and resource allocation.
Operations Research: Used in various industries to optimize processes, reduce costs, and improve efficiency.
Academic Research: Widely used in academic research to develop and test new optimization algorithms and techniques.
Industrial Applications: Applied in industries such as manufacturing, logistics, and telecommunications to solve real-world optimization problems.

By leveraging the strengths of both OR-Tools and SCIP, we can effectively tackle the VRP with time windows and gain valuable insights into the performance and usability of these powerful optimization tools.

Implementation Details

Google OR-Tools Implementation

I modified code from
Vehicle Routing Problem with Time Windows and add capacity constraints to slightly increase complexity and located my code in the file named or-tools.py in the same directory with this file. For mathematician or operations research analyst, this may sound weird. When you write code with OR-tools, you don't formulate mathematical model to code. You just must read the documents and understand how it defines time, distance as cumulative variables and how to set constraints for nodes in a route. Please take a look at the here while reading this explanation.

"""Capacited Vehicles Routing Problem (CVRP) with Time Windows."""

from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp


def create_data_model():
    """Stores the data for the problem."""
    data = {}
    data["time_matrix"] = [
        [0, 6, 9, 8, 7, 3, 6, 2, 3, 2, 6, 6, 4, 4, 5, 9, 7],
        [6, 0, 8, 3, 2, 6, 8, 4, 8, 8, 13, 7, 5, 8, 12, 10, 14],
        [9, 8, 0, 11, 10, 6, 3, 9, 5, 8, 4, 15, 14, 13, 9, 18, 9],
        [8, 3, 11, 0, 1, 7, 10, 6, 10, 10, 14, 6, 7, 9, 14, 6, 16],
        [7, 2, 10, 1, 0, 6, 9, 4, 8, 9, 13, 4, 6, 8, 12, 8, 14],
        [3, 6, 6, 7, 6, 0, 2, 3, 2, 2, 7, 9, 7, 7, 6, 12, 8],
        [6, 8, 3, 10, 9, 2, 0, 6, 2, 5, 4, 12, 10, 10, 6, 15, 5],
        [2, 4, 9, 6, 4, 3, 6, 0, 4, 4, 8, 5, 4, 3, 7, 8, 10],
        [3, 8, 5, 10, 8, 2, 2, 4, 0, 3, 4, 9, 8, 7, 3, 13, 6],
        [2, 8, 8, 10, 9, 2, 5, 4, 3, 0, 4, 6, 5, 4, 3, 9, 5],
        [6, 13, 4, 14, 13, 7, 4, 8, 4, 4, 0, 10, 9, 8, 4, 13, 4],
        [6, 7, 15, 6, 4, 9, 12, 5, 9, 6, 10, 0, 1, 3, 7, 3, 10],
        [4, 5, 14, 7, 6, 7, 10, 4, 8, 5, 9, 1, 0, 2, 6, 4, 8],
        [4, 8, 13, 9, 8, 7, 10, 3, 7, 4, 8, 3, 2, 0, 4, 5, 6],
        [5, 12, 9, 14, 12, 6, 6, 7, 3, 3, 4, 7, 6, 4, 0, 9, 2],
        [9, 10, 18, 6, 8, 12, 15, 8, 13, 9, 13, 3, 4, 5, 9, 0, 9],
        [7, 14, 9, 16, 14, 8, 5, 10, 6, 5, 4, 10, 8, 6, 2, 9, 0],
    ]
    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
        (10, 16),  # 10
        (10, 15),  # 11
        (0, 5),  # 12
        (5, 10),  # 13
        (7, 8),  # 14
        (10, 15),  # 15
        (11, 15),  # 16
    ]
    data["demands"] = [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
    data["vehicle_capacities"] = [15, 15, 15, 15]
    data["num_vehicles"] = 4
    data["depot"] = 0
    return data


def print_solution(data, manager, routing, solution):
    """Prints solution on console."""
    print(f"Objective: {solution.ObjectiveValue()}")
    time_dimension = routing.GetDimensionOrDie("Time")
    total_time = 0
    for vehicle_id in range(data["num_vehicles"]):
        index = routing.Start(vehicle_id)
        plan_output = f"Route for vehicle {vehicle_id}:\n"
        while not routing.IsEnd(index):
            time_var = time_dimension.CumulVar(index)
            plan_output += (
                f"{manager.IndexToNode(index)}"
                f" Time({solution.Min(time_var)},{solution.Max(time_var)})"
                " -> "
            )
            index = solution.Value(routing.NextVar(index))
        time_var = time_dimension.CumulVar(index)
        plan_output += (
            f"{manager.IndexToNode(index)}"
            f" Time({solution.Min(time_var)},{solution.Max(time_var)})\n"
        )
        plan_output += f"Time of the route: {solution.Min(time_var)}min\n"
        print(plan_output)
        total_time += solution.Min(time_var)
    print(f"Total time of all routes: {total_time}min")


def main():
    """Solve the VRP with time windows."""
    # Instantiate the data problem.
    data = create_data_model()

    # 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 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",
    )

    # 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)

    # Print solution on console.
    if solution:
        print_solution(data, manager, routing, solution)


if __name__ == "__main__":
    main()

The code consists of many components, I will explain each component step-by-step.

Routing Index Manager

manager = pywrapcp.RoutingIndexManager(
    len(data["time_matrix"]), data["num_vehicles"], data["depot"]
)

The RoutingIndexManager is responsible for managing the indexing of nodes (locations) and vehicles. It creates a mapping between the problem nodes and internal indices used by the solver. This mapping is crucial for retrieving parameters/values from the data dictionary, which is indexed by node numbers (indices of the numpy array) rather than the internal indices used by Google OR-Tools.

Creating the Routing Model

routing = pywrapcp.RoutingModel(manager)

The RoutingModel is the core of the routing solver. It uses the RoutingIndexManager to manage the nodes and vehicles. When you initialize this object, think of it as declaring the decision variable $x_{ijk}$ , which represents the decision of whether a vehicle $k$ moves from node $i$ to node $j$ .

Registering the Time Callback

def time_callback(from_index, to_index):
    """Returns the travel time between the two nodes."""
    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)

The time_callback function returns the travel time between two nodes. The RegisterTransitCallback method registers this callback with the routing model, allowing it to be used for calculating the travel time between nodes. The transit_callback_index is an identifier for this callback.

Setting the Objective Function

routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

This line sets the objective function as the sum of all travel times for the routes taken by the vehicles. Mathematically, it can be expressed as:

\sum_{i \in N} \sum_{j \in N} \sum_{k \in V} c_{ij} \cdot x_{ijk}

where $c_{ij}$ is the travel time from node $i$ to node $j$ , and $x_{ijk}$ is the decision variable indicating whether vehicle $k$ moves from node $i$ to node $j$ .

Registering the Demand Callback

def demand_callback(from_index):
    """Returns the demand of the node."""
    from_node = manager.IndexToNode(from_index)
    return data["demands"][from_node]

demand_callback_index = routing.RegisterUnaryTransitCallback(demand_callback)

The demand_callback function returns the demand of a node. The RegisterUnaryTransitCallback method registers this callback with the routing model, allowing it to be used for calculating the cumulative demand at each node. The demand_callback_index is an identifier for this callback.

Adding the Capacity Constraint

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

This line adds a capacity constraint to the routing model. The AddDimensionWithVehicleCapacity method sets up a rule to cumulate the number of delivered units at each node. The cumulative variable for each vehicle cannot exceed its capacity limit, specified by data["vehicle_capacities"].

Adding the 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)

This section adds a time windows constraint to the routing model. The AddDimension method sets up a rule to track the cumulative travel time for each vehicle. The parameters are as follows:

transit_callback_index: The callback function for calculating travel time.
30: The maximum allowed waiting time.
30: The maximum travel time per vehicle.
False: Do not force the start cumulative value to zero.
time: The name of the dimension.

The GetDimensionOrDie method retrieves the time dimension, which is used to add time window constraints for each location.

Adding Time Window Constraints for Each Location

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])

This loop adds time window constraints for each location except the depot. The SetRange method sets the allowable time window for each location, ensuring that deliveries are made within the specified time windows.

Adding 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]
    )

This loop adds time window constraints for the start node of each vehicle. The SetRange method ensures that each vehicle starts its route within the specified time window for the depot.

Instantiating Route Start and End Times

for i in range(data["num_vehicles"]):
    routing.AddVariableMinimizedByFinalizer(
        time_dimension.CumulVar(routing.Start(i))
    )
    routing.AddVariableMinimizedByFinalizer(time_dimension.CumulVar(routing.End(i)))

These lines instantiate the start and end times for each vehicle's route, ensuring that the solver produces feasible times for the routes.

Optimizing

# 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)

# Print solution on console.
if solution:
    print_solution(data, manager, routing, solution)

This section sets the search parameters for the solver, specifying the first solution heuristic as PATH_CHEAPEST_ARC. The SolveWithParameters method solves the problem using the specified search parameters. If a solution is found, the print_solution function is called to print the solution on the console.

SCIP Implementations

Node-Based Model

Here is the mathematical formulation, I implemented in this code.

from pyscipopt import Model, quicksum
import networkx as nx

def create_data_model():
    """Stores the data for the problem."""
    data = {}
    # Time matrix representing travel times between locations (consistent with time windows)
    data["time_matrix"] = [
        [0, 6, 9, 8, 7, 3, 6, 2, 3, 2, 6, 6, 4, 4, 5, 9, 7],
        [6, 0, 8, 3, 2, 6, 8, 4, 8, 8, 13, 7, 5, 8, 12, 10, 14],
        [9, 8, 0, 11, 10, 6, 3, 9, 5, 8, 4, 15, 14, 13, 9, 18, 9],
        [8, 3, 11, 0, 1, 7, 10, 6, 10, 10, 14, 6, 7, 9, 14, 6, 16],
        [7, 2, 10, 1, 0, 6, 9, 4, 8, 9, 13, 4, 6, 8, 12, 8, 14],
        [3, 6, 6, 7, 6, 0, 2, 3, 2, 2, 7, 9, 7, 7, 6, 12, 8],
        [6, 8, 3, 10, 9, 2, 0, 6, 2, 5, 4, 12, 10, 10, 6, 15, 5],
        [2, 4, 9, 6, 4, 3, 6, 0, 4, 4, 8, 5, 4, 3, 7, 8, 10],
        [3, 8, 5, 10, 8, 2, 2, 4, 0, 3, 4, 9, 8, 7, 3, 13, 6],
        [2, 8, 8, 10, 9, 2, 5, 4, 3, 0, 4, 6, 5, 4, 3, 9, 5],
        [6, 13, 4, 14, 13, 7, 4, 8, 4, 4, 0, 10, 9, 8, 4, 13, 4],
        [6, 7, 15, 6, 4, 9, 12, 5, 9, 6, 10, 0, 1, 3, 7, 3, 10],
        [4, 5, 14, 7, 6, 7, 10, 4, 8, 5, 9, 1, 0, 2, 6, 4, 8],
        [4, 8, 13, 9, 8, 7, 10, 3, 7, 4, 8, 3, 2, 0, 4, 5, 6],
        [5, 12, 9, 14, 12, 6, 6, 7, 3, 3, 4, 7, 6, 4, 0, 9, 2],
        [9, 10, 18, 6, 8, 12, 15, 8, 13, 9, 13, 3, 4, 5, 9, 0, 9],
        [7, 14, 9, 16, 14, 8, 5, 10, 6, 5, 4, 10, 8, 6, 2, 9, 0],
    ]
    data["time_windows"] = [
        (0, 5),     # 0 - 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
        (10, 16),   # 10
        (10, 15),   # 11
        (0, 5),     # 12
        (5, 10),    # 13
        (7, 8),     # 14
        (10, 15),   # 15
        (11, 15),   # 16
    ]
    data["demands"] = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]
    data["vehicle_capacities"] = [15, 15, 15, 15]  # Capacities for each vehicle
    data["num_vehicles"] = 4
    data["depot"] = 0
    return data

def main():
    """Solve the VRPTW problem using PySCIPOpt with variable depot departure times."""
    # Retrieve data
    data = create_data_model()
    time_matrix = data["time_matrix"]
    time_windows = data["time_windows"]
    demands = data["demands"]
    vehicle_capacities = data["vehicle_capacities"]
    num_vehicles = data["num_vehicles"]
    depot = data["depot"]

    n = len(time_matrix)  # Number of nodes
    N = range(n)  # Set of nodes
    V = [i for i in N if i != depot]  # Customers excluding depot
    K = range(num_vehicles)  # Set of vehicles

    # Calculate M for Big-M constraints (maximum possible time)
    max_time = max(tw[1] for tw in time_windows) + max(
        time_matrix[i][j] for i in N for j in N
    )

    # Initialize the model
    model = Model("VRPTW with Variable Depot Departure Times")

    # Variables
    x = {}  # Binary variables: x[i,j,k] == 1 if vehicle k travels from i to j
    t = {}  # Arrival time at node i by any vehicle
    t_start = {}  # Departure time from depot for vehicle k
    t_end = {}    # Return time to depot for vehicle k
    w = {}  # Waiting time at node i by vehicle k
    y = {}  # Binary variables: y[i,k] == 1 if customer i is visited by vehicle k

    # Decision variables
    for k in K:
        capacity_k = vehicle_capacities[k]
        # Departure time from depot
        t_start[k] = model.addVar(lb=time_windows[depot][0], ub=time_windows[depot][1], vtype="I", name=f"t_start_{k}")
        # Return time to depot
        t_end[k] = model.addVar(lb=time_windows[depot][0], ub=None, vtype="I", name=f"t_end_{k}")
        for i in N:
            if i != depot:
                y[i, k] = model.addVar(vtype="B", name=f"y_{i}_{k}")
            for j in N:
                if i != j:
                    x[i, j, k] = model.addVar(vtype="B", name=f"x_{i}_{j}_{k}")

    for i in N:
        if i != depot:
            t[i] = model.addVar(lb=time_windows[i][0], ub=time_windows[i][1], vtype="I", name=f"t_{i}")

    for i in V:
        for k in K:
            w[i, k] = model.addVar(lb=0, ub=max_time, vtype="I", name=f"w_{i}_{k}")

    # Objective Function: Minimize total operation time of all vehicles
    model.setObjective(
        quicksum(t_end[k] - t_start[k] for k in K),
        "minimize",
    )

    # Constraints

    # 1. Each customer is visited exactly once
    for i in V:
        model.addCons(
            quicksum(y[i, k] for k in K) == 1,
            name=f"VisitOnce_{i}"
        )

    # 2. For each vehicle, ensure it departs from and returns to the depot correctly
    for k in K:
        model.addCons(
            quicksum(x[depot, j, k] for j in N if j != depot) == 1,
            name=f"DepartDepot_{k}"
        )
        model.addCons(
            quicksum(x[i, depot, k] for i in N if i != depot) == 1,
            name=f"ReturnDepot_{k}"
        )

    # 3. Link x and y variables
    for k in K:
        for i in V:
            model.addCons(
                quicksum(x[i, j, k] for j in N if j != i) == y[i, k],
                name=f"LinkXY_{i}_{k}"
            )
            model.addCons(
                quicksum(x[j, i, k] for j in N if j != i) == y[i, k],
                name=f"LinkXY_Inflow_{i}_{k}"
            )

    # 4. Flow conservation constraints for nodes (excluding depot)
    for k in K:
        for i in V:
            inflow = quicksum(x[j, i, k] for j in N if j != i)
            outflow = quicksum(x[i, j, k] for j in N if j != i)
            model.addCons(
                inflow - outflow == 0,
                name=f"FlowConservation_{i}_{k}"
            )

    # 5. Capacity constraints for each vehicle
    for k in K:
        capacity_k = vehicle_capacities[k]
        model.addCons(
            quicksum(demands[i] * y[i, k] for i in V) <= capacity_k,
            name=f"VehicleCapacity_{k}"
        )

    # 6. Time window constraints and time propagation constraints
    for k in K:
        # Departure time from depot is within time window
        model.addCons(
            t_start[k] >= time_windows[depot][0],
            name=f"DepotTimeWindowStart_{k}"
        )
        model.addCons(
            t_start[k] <= time_windows[depot][1],
            name=f"DepotTimeWindowEnd_{k}"
        )

        # Time propagation constraints
        for i in N:
            for j in N:
                if i != j:
                    if i == depot and j != depot:
                        # From depot to first customer
                        model.addCons(
                            t[j] >= t_start[k] + time_matrix[i][j] - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagationDepot_{i}_{j}_{k}"
                        )
                    elif i != depot and j != depot:
                        # Between customers
                        model.addCons(
                            t[j] >= t[i] + w[i, k] + time_matrix[i][j] - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagation_{i}_{j}_{k}"
                        )
                    elif i != depot and j == depot:
                        # Return to depot
                        model.addCons(
                            t_end[k] >= t[i] + w[i, k] + time_matrix[i][j] - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagationReturnDepot_{i}_{j}_{k}"
                        )

        # Waiting time constraints
        for i in V:
            # Waiting time is non-negative
            model.addCons(
                w[i, k] >= 0,
                name=f"WaitingTimeNonNegative_{i}_{k}"
            )
            # If customer i is not visited by vehicle k, waiting time is zero
            model.addCons(
                w[i, k] <= max_time * y[i, k],
                name=f"WaitingTimeZeroIfNotVisited_{i}_{k}"
            )

    # Time windows for customer arrivals
    for i in V:
        model.addCons(
            t[i] >= time_windows[i][0],
            name=f"TimeWindowStart_{i}"
        )
        model.addCons(
            t[i] <= time_windows[i][1],
            name=f"TimeWindowEnd_{i}"
        )

    # Waiting time zero if not visited (already added)

    # Solve the model
    model.optimize()

    # Retrieve and print the solution
    status = model.getStatus()
    if status == "optimal" or status == "bestsollimit":
        print_solution(data, x, t, t_start, t_end, w, K, N, depot, model)
    else:
        print("No optimal solution found.")
        print(f"Solver status: {status}")

def print_solution(data, x, t, t_start, t_end, w, K, N, depot, model):
    """Prints solution in the specified format."""
    solution = model.getBestSol()
    total_time = 0

    print(f"Objective: {model.getSolObjVal(solution):.1f}")
    for k in K:
        vehicle_id = k
        # Build the route for vehicle k
        edges = [(i, j) for i in N for j in N if i != j and model.getSolVal(solution, x.get((i, j, k), 0)) > 0.5]
        if not edges:
            continue  # Skip if vehicle k is not used
        G = nx.DiGraph()
        G.add_edges_from(edges)

        if depot not in G.nodes:
            continue  # Skip if vehicle k is not used

        plan_output = f"Route for vehicle {vehicle_id}:\n"

        index = depot
        time = model.getSolVal(solution, t_start[k])
        time_str = f"Time({time:.1f},{time:.1f})"
        plan_output += f"{index} {time_str} -> "
        total_route_time = time

        while True:
            successors = list(G.successors(index))
            if not successors or successors[0] == depot:
                break
            next_node = successors[0]
            arrival_time = model.getSolVal(solution, t[next_node])
            arrival_time_str = f"Time({arrival_time:.1f},{arrival_time:.1f})"
            plan_output += f"{next_node} {arrival_time_str} -> "
            total_route_time = arrival_time
            index = next_node

        # Return to depot
        return_time = model.getSolVal(solution, t_end[k])
        arrival_time_str = f"Time({return_time:.1f},{return_time:.1f})"
        plan_output += f"{depot} {arrival_time_str}\n"
        route_duration = return_time - model.getSolVal(solution, t_start[k])
        plan_output += f"Time of the route: {route_duration:.1f}min\n"
        print(plan_output)
        total_time += route_duration

    print(f"Total time of all routes: {total_time:.1f}min")

if __name__ == "__main__":
    main()

Sets and Indices

$N$ : Set of all nodes (including depot)
$V$ : Set of customer nodes (excluding depot)
$K$ : Set of vehicles
$i, j \in N$ : Indices for nodes
$k \in K$ : Index for vehicles

Parameters

$c_{ij}$ : Travel time from node $i$ to node $j$
$[a_i, b_i]$ : Time window for node $i$
$d_i$ : Demand at node $i$
$Q_k$ : Capacity of vehicle $k$
$M$ : A large constant (Big-M)

Decision Variables

$x_{ijk}$ : Binary variable, 1 if vehicle $k$ travels from node $i$ to node $j$ , 0 otherwise
$t_i$ : Arrival time at node $i$
$t_{start_k}$ : Departure time from depot for vehicle $k$
$t_{end_k}$ : Return time to depot for vehicle $k$
$w_{ik}$ : Waiting time at node $i$ by vehicle $k$
$y_{ik}$ : Binary variable, 1 if customer $i$ is visited by vehicle $k$ , 0 otherwise

Objective Function

Minimize the total operation time of all vehicles:

\text{Minimize} \quad \sum_{k \in K} (t_{end_k} - t_{start_k})

Constraints

1. Each customer is visited exactly once

\sum_{k \in K} y_{ik} = 1 \quad \forall i \in V

2. Each vehicle departs from and returns to the depot exactly once

\sum_{j \in N, j \neq i} x_{0jk} = 1 \quad \forall k \in K

\sum_{i \in N, i \neq j} x_{i0k} = 1 \quad \forall k \in K

3. Link $x$ and $y$ variables

\sum_{j \in N, j \neq i} x_{ijk} = y_{ik} \quad \forall i \in V, \forall k \in K

\sum_{j \in N, j \neq i} x_{jik} = y_{ik} \quad \forall i \in V, \forall k \in K

4. Flow conservation constraints for nodes (excluding depot)

\sum_{j \in N, j \neq i} x_{jik} - \sum_{j \in N, j \neq i} x_{ijk} = 0 \quad \forall i \in V, \forall k \in K

5. Capacity constraints for each vehicle

\sum_{i \in V} d_i y_{ik} \leq Q_k \quad \forall k \in K

6. Time window constraints and time propagation constraints

Departure time from depot is within time window

a_0 \leq t_{\text{start}_k} \leq b_0 \quad \forall k \in K

Time propagation constraints

From depot to first customer:

$t_j \geq t_{start_k} + c_{0j} - M(1 - x_{0jk}) \quad \forall j \in V, \forall k \in K$
Between customers:

$t_j \geq t_i + w_{ik} + c_{ij} - M(1 - x_{ijk}) \quad \forall i, j \in V, i \neq j, \forall k \in K$
Return to depot:

$t_{end_k} \geq t_i + w_{ik} + c_{i0} - M(1 - x_{i0k}) \quad \forall i \in V, \forall k \in K$

Flow-based Model

Here is the mathematical formulation, I implemented in this code.

from pyscipopt import Model, quicksum
import networkx as nx

def create_data_model():
    """Stores the data for the problem."""
    data = {}
    # Time matrix representing travel times between locations (consistent with time windows)
    data["time_matrix"] = [
        [0, 6, 9, 8, 7, 3, 6, 2, 3, 2, 6, 6, 4, 4, 5, 9, 7],
        [6, 0, 8, 3, 2, 6, 8, 4, 8, 8, 13, 7, 5, 8, 12, 10, 14],
        [9, 8, 0, 11, 10, 6, 3, 9, 5, 8, 4, 15, 14, 13, 9, 18, 9],
        [8, 3, 11, 0, 1, 7, 10, 6, 10, 10, 14, 6, 7, 9, 14, 6, 16],
        [7, 2, 10, 1, 0, 6, 9, 4, 8, 9, 13, 4, 6, 8, 12, 8, 14],
        [3, 6, 6, 7, 6, 0, 2, 3, 2, 2, 7, 9, 7, 7, 6, 12, 8],
        [6, 8, 3, 10, 9, 2, 0, 6, 2, 5, 4, 12, 10, 10, 6, 15, 5],
        [2, 4, 9, 6, 4, 3, 6, 0, 4, 4, 8, 5, 4, 3, 7, 8, 10],
        [3, 8, 5, 10, 8, 2, 2, 4, 0, 3, 4, 9, 8, 7, 3, 13, 6],
        [2, 8, 8, 10, 9, 2, 5, 4, 3, 0, 4, 6, 5, 4, 3, 9, 5],
        [6, 13, 4, 14, 13, 7, 4, 8, 4, 4, 0, 10, 9, 8, 4, 13, 4],
        [6, 7, 15, 6, 4, 9, 12, 5, 9, 6, 10, 0, 1, 3, 7, 3, 10],
        [4, 5, 14, 7, 6, 7, 10, 4, 8, 5, 9, 1, 0, 2, 6, 4, 8],
        [4, 8, 13, 9, 8, 7, 10, 3, 7, 4, 8, 3, 2, 0, 4, 5, 6],
        [5, 12, 9, 14, 12, 6, 6, 7, 3, 3, 4, 7, 6, 4, 0, 9, 2],
        [9, 10, 18, 6, 8, 12, 15, 8, 13, 9, 13, 3, 4, 5, 9, 0, 9],
        [7, 14, 9, 16, 14, 8, 5, 10, 6, 5, 4, 10, 8, 6, 2, 9, 0],
    ]
    data["time_windows"] = [
        (0, 5),     # 0 - 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
        (10, 16),   # 10
        (10, 15),   # 11
        (0, 5),     # 12
        (5, 10),    # 13
        (7, 8),     # 14
        (10, 15),   # 15
        (11, 15),   # 16
    ]
    data["demands"] = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]
    data["vehicle_capacities"] = [15, 15, 15, 15]  # Capacities for each vehicle
    data["num_vehicles"] = 4
    data["depot"] = 0
    return data

def main():
    """Solve the VRPTW problem using PySCIPOpt with a flow-based model."""
    # Retrieve data
    data = create_data_model()
    time_matrix = data["time_matrix"]
    time_windows = data["time_windows"]
    demands = data["demands"]
    vehicle_capacities = data["vehicle_capacities"]
    num_vehicles = data["num_vehicles"]
    depot = data["depot"]

    n = len(time_matrix)  # Number of nodes
    N = range(n)  # Set of nodes
    V = [i for i in N if i != depot]  # Customers excluding depot
    K = range(num_vehicles)  # Set of vehicles

    # Calculate M for Big-M constraints (maximum possible time)
    max_time = max(tw[1] for tw in time_windows) + max(
        time_matrix[i][j] for i in N for j in N
    )

    # Initialize the model
    model = Model("VRPTW with Flow-Based Model")

    # Variables
    x = {}  # Binary variables: x[i,j,k] == 1 if vehicle k travels from i to j
    f = {}  # Continuous variables: flow of demand on arc (i,j) by vehicle k
    t = {}  # Arrival time at node i by any vehicle
    t_start = {}  # Departure time from depot for vehicle k
    t_end = {}    # Return time to depot for vehicle k
    w = {}  # Waiting time at node i by vehicle k
    y = {}  # Binary variables: y[i,k] == 1 if customer i is visited by vehicle k

    # Decision variables
    for k in K:
        capacity_k = vehicle_capacities[k]
        # Departure time from depot
        t_start[k] = model.addVar(lb=time_windows[depot][0], ub=time_windows[depot][1], vtype="I", name=f"t_start_{k}")
        # Return time to depot
        t_end[k] = model.addVar(lb=time_windows[depot][0], ub=None, vtype="I", name=f"t_end_{k}")
        for i in N:
            if i != depot:
                y[i, k] = model.addVar(vtype="B", name=f"y_{i}_{k}")
            for j in N:
                if i != j:
                    x[i, j, k] = model.addVar(vtype="B", name=f"x_{i}_{j}_{k}")
                    # Flow variable
                    f[i, j, k] = model.addVar(lb=0, vtype="I", name=f"f_{i}_{j}_{k}")

    for i in N:
        if i != depot:
            t[i] = model.addVar(lb=time_windows[i][0], ub=time_windows[i][1], vtype="I", name=f"t_{i}")

    for i in V:
        for k in K:
            w[i, k] = model.addVar(lb=0, ub=max_time, vtype="I", name=f"w_{i}_{k}")

    # Objective Function: Minimize total operation time of all vehicles
    model.setObjective(
        quicksum(t_end[k] - t_start[k] for k in K),
        "minimize",
    )

    # Constraints

    # 1. Each customer is visited exactly once
    for i in V:
        model.addCons(
            quicksum(y[i, k] for k in K) == 1,
            name=f"VisitOnce_{i}"
        )

    # 2. For each vehicle, ensure it departs from and returns to the depot correctly
    for k in K:
        model.addCons(
            quicksum(x[depot, j, k] for j in N if j != depot) == 1,
            name=f"DepartDepot_{k}"
        )
        model.addCons(
            quicksum(x[i, depot, k] for i in N if i != depot) == 1,
            name=f"ReturnDepot_{k}"
        )

    # 3. Link x and y variables
    for k in K:
        for i in V:
            model.addCons(
                quicksum(x[i, j, k] for j in N if j != i) == y[i, k],
                name=f"LinkXY_{i}_{k}"
            )
            model.addCons(
                quicksum(x[j, i, k] for j in N if j != i) == y[i, k],
                name=f"LinkXY_Inflow_{i}_{k}"
            )

    # 4. Flow conservation constraints (using flow variables)
    for k in K:
        capacity_k = vehicle_capacities[k]
        # Capacity constraints on arcs
        for i in N:
            for j in N:
                if i != j:
                    model.addCons(
                        f[i, j, k] <= capacity_k * x[i, j, k],
                        name=f"CapacityArc_{i}_{j}_{k}"
                    )
        # Flow conservation at depot
        inflow_f = quicksum(f[j, depot, k] for j in N if j != depot)
        outflow_f = quicksum(f[depot, j, k] for j in N if j != depot)
        total_demand_k = quicksum(demands[i] * y[i, k] for i in V)
        model.addCons(
            inflow_f - outflow_f == -total_demand_k,
            name=f"FlowConservation_Depot_{k}"
        )
        # Flow conservation at customers
        for i in V:
            inflow_f = quicksum(f[j, i, k] for j in N if j != i)
            outflow_f = quicksum(f[i, j, k] for j in N if j != i)
            model.addCons(
                inflow_f - outflow_f == demands[i] * y[i, k],
                name=f"FlowConservation_{i}_{k}"
            )

    # 5. Time window constraints and time propagation constraints
    for k in K:
        # Departure time from depot is within time window
        model.addCons(
            t_start[k] >= time_windows[depot][0],
            name=f"DepotTimeWindowStart_{k}"
        )
        model.addCons(
            t_start[k] <= time_windows[depot][1],
            name=f"DepotTimeWindowEnd_{k}"
        )

        # Time propagation constraints
        for i in N:
            for j in N:
                if i != j:
                    if i == depot and j != depot:
                        # From depot to first customer
                        model.addCons(
                            t[j] >= t_start[k] + time_matrix[i][j] - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagationDepot_{i}_{j}_{k}"
                        )
                    elif i != depot and j != depot:
                        # Between customers
                        model.addCons(
                            t[j] >= t[i] + w[i, k] + time_matrix[i][j] - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagation_{i}_{j}_{k}"
                        )
                    elif i != depot and j == depot:
                        # Return to depot
                        model.addCons(
                            t_end[k] >= t[i] + w[i, k] + time_matrix[i][j] - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagationReturnDepot_{i}_{j}_{k}"
                        )

        # Waiting time constraints
        for i in V:
            # Waiting time is non-negative
            model.addCons(w[i, k] >= 0,
                name=f"WaitingTimeNonNegative_{i}_{k}"
            )
            # If customer i is not visited by vehicle k, waiting time is zero
            model.addCons(
                w[i, k] <= max_time * y[i, k],
                name=f"WaitingTimeZeroIfNotVisited_{i}_{k}"
            )

    # Time windows for customer arrivals
    for i in V:
        model.addCons(
            t[i] >= time_windows[i][0],
            name=f"TimeWindowStart_{i}"
        )
        model.addCons(
            t[i] <= time_windows[i][1],
            name=f"TimeWindowEnd_{i}"
        )

    # Solve the model
    model.optimize()

    # Retrieve and print the solution
    status = model.getStatus()
    if status == "optimal" or status == "bestsollimit":
        print_solution(data, x, t, t_start, t_end, w, K, N, depot, model)
    else:
        print("No optimal solution found.")
        print(f"Solver status: {status}")

def print_solution(data, x, t, t_start, t_end, w, K, N, depot, model):
    """Prints solution in the specified format."""
    solution = model.getBestSol()
    total_time = 0

    print(f"Objective: {model.getSolObjVal(solution):.1f}")
    for k in K:
        vehicle_id = k
        # Build the route for vehicle k
        edges = [(i, j) for i in N for j in N if i != j and model.getSolVal(solution, x.get((i, j, k), 0)) > 0.5]
        if not edges:
            continue  # Skip if vehicle k is not used
        G = nx.DiGraph()
        G.add_edges_from(edges)

        if depot not in G.nodes:
            continue  # Skip if vehicle k is not used

        plan_output = f"Route for vehicle {vehicle_id}:\n"

        index = depot
        time = model.getSolVal(solution, t_start[k])
        time_str = f"Time({time:.1f},{time:.1f})"
        plan_output += f"{index} {time_str} -> "
        total_route_time = time

        while True:
            successors = list(G.successors(index))
            if not successors or successors[0] == depot:
                break
            next_node = successors[0]
            arrival_time = model.getSolVal(solution, t[next_node])
            arrival_time_str = f"Time({arrival_time:.1f},{arrival_time:.1f})"
            plan_output += f"{next_node} {arrival_time_str} -> "
            total_route_time = arrival_time
            index = next_node

        # Return to depot
        return_time = model.getSolVal(solution, t_end[k])
        arrival_time_str = f"Time({return_time:.1f},{return_time:.1f})"
        plan_output += f"{depot} {arrival_time_str}\n"
        route_duration = return_time - model.getSolVal(solution, t_start[k])
        plan_output += f"Time of the route: {route_duration:.1f}min\n"
        print(plan_output)
        total_time += route_duration

    print(f"Total time of all routes: {total_time:.1f}min")

if __name__ == "__main__":
    main()

Sets and Indices

$N$ : Set of all nodes (including depot)
$V$ : Set of customer nodes (excluding depot)
$K$ : Set of vehicles
$i, j \in N$ : Indices for nodes
$k \in K$ : Index for vehicles

Parameters

$c_{ij}$ : Travel time from node $i$ to node $j$
$[a_i, b_i]$ : Time window for node $i$
$d_i$ : Demand at node $i$
$Q_k$ : Capacity of vehicle $k$
$M$ : A large constant (Big-M)

Decision Variables

$x_{ijk}$ : Binary variable, 1 if vehicle $k$ travels from node $i$ to node $j$ , 0 otherwise
$f_{ijk}$ : Continuous variable, flow of demand on arc $(i, j)$ by vehicle $k$
$t_i$ : Arrival time at node $i$
$t_{start_k}$ : Departure time from depot for vehicle $k$
$t_{end_k}$ : Return time to depot for vehicle $k$
$w_{ik}$ : Waiting time at node $i$ by vehicle $k$
$y_{ik}$ : Binary variable, 1 if customer $i$ is visited by vehicle $k$ , 0 otherwise

Objective Function

Minimize the total operation time of all vehicles:

\text{Minimize} \quad \sum_{k \in K} (t_{end_k} - t_{start_k})

Constraints

1. Each customer is visited exactly once

\sum_{k \in K} y_{ik} = 1 \quad \forall i \in V

2. Each vehicle departs from and returns to the depot exactly once

\sum_{j \in N, j \neq i} x_{0jk} = 1 \quad \forall k \in K

\sum_{i \in N, i \neq j} x_{i0k} = 1 \quad \forall k \in K

3. Link $x$ and $y$ variables

\sum_{j \in N, j \neq i} x_{ijk} = y_{ik} \quad \forall i \in V, \forall k \in K

\sum_{j \in N, j \neq i} x_{jik} = y_{ik} \quad \forall i \in V, \forall k \in K

4. Flow conservation constraints (using flow variables)

Capacity constraints on arcs:

$f_{ijk} \leq Q_k \cdot x_{ijk} \quad \forall i, j \in N, i \neq j, \forall k \in K$
Flow conservation at depot:

$\sum_{j \in N, j \neq 0} f_{j0k} - \sum_{j \in N, j \neq 0} f_{0jk} = -\sum_{i \in V} d_i \cdot y_{ik} \quad \forall k \in K$
Flow conservation at customers:

$\sum_{j \in N, j \neq i} f_{jik} - \sum_{j \in N, j \neq i} f_{ijk} = d_i \cdot y_{ik} \quad \forall i \in V, \forall k \in K$

5. Time window constraints and time propagation constraints

Departure time from depot is within time window:

$a_0 \leq t_{start_k} \leq b_0 \quad \forall k \in K$
Time propagation constraints:
- From depot to first customer: $t_j \geq t_{start_k} + c_{0j} - M(1 - x_{0jk}) \quad \forall j \in V, \forall k \in K$
- Between customers: $t_j \geq t_i + w_{ik} + c_{ij} - M(1 - x_{ijk}) \quad \forall i, j \in V, i \neq j, \forall k \in K$
- Return to depot: $t_{end_k} \geq t_i + w_{ik} + c_{i0} - M(1 - x_{i0k}) \quad \forall i \in V, \forall k \in K$
Waiting time constraints:
- Waiting time is non-negative: $w_{ik} \geq 0 \quad \forall i \in V, \forall k \in K$
- If customer $i$ is not visited by vehicle $k$ , waiting time is zero: $w_{ik} \leq M \cdot y_{ik} \quad \forall i \in V, \forall k \in K$

6. Time windows for customer arrivals

a_i \leq t_i \leq b_i \quad \forall i \in V

MTZ Model

Here is the mathematical formulation, I implemented in this code.

from pyscipopt import Model, quicksum
import networkx as nx

def create_data_model():
    """Stores the data for the problem."""
    data = {}
    # Time matrix representing travel times between locations (consistent with time windows)
    data["time_matrix"] = [
        [0, 6, 9, 8, 7, 3, 6, 2, 3, 2, 6, 6, 4, 4, 5, 9, 7],
        [6, 0, 8, 3, 2, 6, 8, 4, 8, 8, 13, 7, 5, 8, 12, 10, 14],
        [9, 8, 0, 11, 10, 6, 3, 9, 5, 8, 4, 15, 14, 13, 9, 18, 9],
        [8, 3, 11, 0, 1, 7, 10, 6, 10, 10, 14, 6, 7, 9, 14, 6, 16],
        [7, 2, 10, 1, 0, 6, 9, 4, 8, 9, 13, 4, 6, 8, 12, 8, 14],
        [3, 6, 6, 7, 6, 0, 2, 3, 2, 2, 7, 9, 7, 7, 6, 12, 8],
        [6, 8, 3, 10, 9, 2, 0, 6, 2, 5, 4, 12, 10, 10, 6, 15, 5],
        [2, 4, 9, 6, 4, 3, 6, 0, 4, 4, 8, 5, 4, 3, 7, 8, 10],
        [3, 8, 5, 10, 8, 2, 2, 4, 0, 3, 4, 9, 8, 7, 3, 13, 6],
        [2, 8, 8, 10, 9, 2, 5, 4, 3, 0, 4, 6, 5, 4, 3, 9, 5],
        [6, 13, 4, 14, 13, 7, 4, 8, 4, 4, 0, 10, 9, 8, 4, 13, 4],
        [6, 7, 15, 6, 4, 9, 12, 5, 9, 6, 10, 0, 1, 3, 7, 3, 10],
        [4, 5, 14, 7, 6, 7, 10, 4, 8, 5, 9, 1, 0, 2, 6, 4, 8],
        [4, 8, 13, 9, 8, 7, 10, 3, 7, 4, 8, 3, 2, 0, 4, 5, 6],
        [5, 12, 9, 14, 12, 6, 6, 7, 3, 3, 4, 7, 6, 4, 0, 9, 2],
        [9, 10, 18, 6, 8, 12, 15, 8, 13, 9, 13, 3, 4, 5, 9, 0, 9],
        [7, 14, 9, 16, 14, 8, 5, 10, 6, 5, 4, 10, 8, 6, 2, 9, 0],
    ]
    data["time_windows"] = [
        (0, 5),     # 0 - 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
        (10, 16),   # 10
        (10, 15),   # 11
        (0, 5),     # 12
        (5, 10),    # 13
        (7, 8),     # 14
        (10, 15),   # 15
        (11, 15),   # 16
    ]
    data["demands"] = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]
    data["vehicle_capacities"] = [15, 15, 15, 15]  # Capacities for each vehicle
    data["num_vehicles"] = 4
    data["depot"] = 0
    return data

def main():
    """Solve the VRPTW problem using PySCIPOpt with MTZ constraints."""
    # Retrieve data
    data = create_data_model()
    time_matrix = data["time_matrix"]
    time_windows = data["time_windows"]
    demands = data["demands"]
    vehicle_capacities = data["vehicle_capacities"]
    num_vehicles = data["num_vehicles"]
    depot = data["depot"]

    n = len(time_matrix)  # Number of nodes
    N = range(n)  # Set of nodes
    V = [i for i in N if i != depot]  # Customers excluding depot
    K = range(num_vehicles)  # Set of vehicles

    # Calculate M for Big-M constraints (maximum possible time)
    max_time = max(tw[1] for tw in time_windows) + max(
        time_matrix[i][j] for i in N for j in N
    )

    # Initialize the model
    model = Model("VRPTW with MTZ Constraints")

    # Variables
    x = {}  # Binary variables: x[i,j,k] == 1 if vehicle k travels from i to j
    t = {}  # Arrival time at node i by vehicle k
    t_start = {}  # Departure time from depot for vehicle k
    t_end = {}    # Return time to depot for vehicle k
    y = {}  # Binary variables: y[i,k] == 1 if vehicle k visits node i
    u = {}  # MTZ variables for subtour elimination

    # Decision variables
    for k in K:
        capacity_k = vehicle_capacities[k]
        # Departure time from depot
        t_start[k] = model.addVar(lb=time_windows[depot][0], ub=time_windows[depot][1], vtype="I", name=f"t_start_{k}")
        # Return time to depot
        t_end[k] = model.addVar(lb=time_windows[depot][0], ub=None, vtype="I", name=f"t_end_{k}")
        u[depot, k] = model.addVar(lb=0, ub=n, vtype="I", name=f"u_{depot}_{k}")  # u at depot is 0
        model.addCons(u[depot, k] == 0, name=f"MTZ_depot_{k}")
        for i in N:
            if i != depot:
                y[i, k] = model.addVar(vtype="B", name=f"y_{i}_{k}")
                u[i, k] = model.addVar(lb=1, ub=n, vtype="I", name=f"u_{i}_{k}")
                t[i, k] = model.addVar(lb=time_windows[i][0], ub=time_windows[i][1], vtype="I", name=f"t_{i}_{k}")
            for j in N:
                if i != j:
                    x[i, j, k] = model.addVar(vtype="B", name=f"x_{i}_{j}_{k}")

    # Objective Function: Minimize total operation time of all vehicles
    model.setObjective(
        quicksum(t_end[k] - t_start[k] for k in K),
        "minimize",
    )

    # Constraints

    # 1. Each customer is visited exactly once
    for i in V:
        model.addCons(
            quicksum(y[i, k] for k in K) == 1,
            name=f"VisitOnce_{i}"
        )

    # 2. For each vehicle, ensure it departs from and returns to the depot correctly
    for k in K:
        model.addCons(
            quicksum(x[depot, j, k] for j in N if j != depot) == 1,
            name=f"DepartDepot_{k}"
        )
        model.addCons(
            quicksum(x[i, depot, k] for i in N if i != depot) == 1,
            name=f"ReturnDepot_{k}"
        )

    # 3. Flow conservation constraints for nodes (excluding depot)
    for k in K:
        for i in V:
            model.addCons(
                quicksum(x[i, j, k] for j in N if j != i) - quicksum(x[j, i, k] for j in N if j != i) == 0,
                name=f"FlowConservation_{i}_{k}"
            )

    # 4. Link x and y variables
    for k in K:
        for i in V:
            model.addCons(
                quicksum(x[i, j, k] for j in N if j != i) == y[i, k],
                name=f"LinkXY_{i}_{k}"
            )

    # 5. Capacity constraints for each vehicle
    for k in K:
        capacity_k = vehicle_capacities[k]
        model.addCons(
            quicksum(demands[i] * y[i, k] for i in V) <= capacity_k,
            name=f"VehicleCapacity_{k}"
        )

    # 6. Time window constraints and time propagation constraints
    for k in K:
        # Departure time from depot is within time window
        model.addCons(
            t_start[k] >= time_windows[depot][0],
            name=f"DepotTimeWindowStart_{k}"
        )
        model.addCons(
            t_start[k] <= time_windows[depot][1],
            name=f"DepotTimeWindowEnd_{k}"
        )

        # Time propagation constraints
        for i in N:
            if i != depot:
                # Arrival time at node i for vehicle k
                model.addCons(
                    t[i, k] >= time_windows[i][0],
                    name=f"TimeWindowStart_{i}_{k}"
                )
                model.addCons(
                    t[i, k] <= time_windows[i][1],
                    name=f"TimeWindowEnd_{i}_{k}"
                )
            for j in N:
                if i != j:
                    travel_time = time_matrix[i][j]
                    if i == depot:
                        # From depot to first customer
                        model.addCons(
                            t[j, k] >= t_start[k] + travel_time - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagationDepot_{i}_{j}_{k}"
                        )
                    elif j == depot:
                        # From last customer to depot
                        model.addCons(
                            t_end[k] >= t[i, k] + travel_time - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagationReturnDepot_{i}_{j}_{k}"
                        )
                    else:
                        # Between customers
                        model.addCons(
                            t[j, k] >= t[i, k] + travel_time - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagation_{i}_{j}_{k}"
                        )

    # 7. Subtour elimination (MTZ) constraints
    for k in K:
        for i in N:
            if i != depot:
                model.addCons(
                    u[i, k] >= 1,
                    name=f"MTZ_LB_{i}_{k}"
                )
                model.addCons(
                    u[i, k] <= n,
                    name=f"MTZ_UB_{i}_{k}"
                )
        for i in N:
            for j in N:
                if i != j and i != depot and j != depot:
                    model.addCons(
                        u[i, k] - u[j, k] + n * x[i, j, k] <= n - 1,
                        name=f"MTZ_{i}_{j}_{k}"
                    )

    # 8. Ensure that if a customer is not visited by vehicle k, its time variables are inactive
    for k in K:
        for i in V:
            model.addCons(
                y[i, k] <= quicksum(x[j, i, k] for j in N if j != i),
                name=f"VisitIndicator_{i}_{k}"
            )

    # Solve the model
    model.optimize()

    # Retrieve and print the solution
    status = model.getStatus()
    if status == "optimal" or status == "bestsollimit":
        print_solution(data, x, t_start, t_end, t, y, K, N, depot, model)
    else:
        print("No optimal solution found.")
        print(f"Solver status: {status}")

def print_solution(data, x, t_start, t_end, t, y, K, N, depot, model):
    """Prints solution in the specified format."""
    solution = model.getBestSol()
    total_time = 0

    print(f"Objective: {model.getSolObjVal(solution):.1f}")
    for k in K:
        vehicle_id = k
        # Build the route for vehicle k
        edges = [(i, j) for i in N for j in N if i != j and model.getSolVal(solution, x.get((i, j, k), 0)) > 0.5]
        if not edges:
            continue  # Skip if vehicle k is not used
        G = nx.DiGraph()
        G.add_edges_from(edges)

        if depot not in G.nodes:
            continue  # Skip if vehicle k is not used

        plan_output = f"Route for vehicle {vehicle_id}:\n"

        route = [depot]
        time_log = [f"Time({model.getSolVal(solution, t_start[k]):.1f})"]
        next_node = depot

        while True:
            successors = list(G.successors(next_node))
            if not successors:
                break
            next_node = successors[0]
            if next_node == depot:
                break
            route.append(next_node)
            arrival_time = model.getSolVal(solution, t[next_node, k])
            time_log.append(f"Time({arrival_time:.1f})")

        # Append depot at the end
        route.append(depot)
        time_log.append(f"Time({model.getSolVal(solution, t_end[k]):.1f})")

        # Print the route and times
        for idx, node in enumerate(route):
            plan_output += f"{node} {time_log[idx]} -> "
        plan_output = plan_output.rstrip(" -> ") + "\n"
        route_duration = model.getSolVal(solution, t_end[k]) - model.getSolVal(solution, t_start[k])
        plan_output += f"Time of the route: {route_duration:.1f}min\n"
        print(plan_output)
        total_time += route_duration

    print(f"Total time of all routes: {total_time:.1f}min")

if __name__ == "__main__":
    main()

Sets and Indices

$N$ : Set of all nodes (including depot)
$V$ : Set of customer nodes (excluding depot)
$K$ : Set of vehicles
$i, j \in N$ : Indices for nodes
$k \in K$ : Index for vehicles

Parameters

$c_{ij}$ : Travel time from node $i$ to node $j$
$[a_i, b_i]$ : Time window for node $i$
$d_i$ : Demand at node $i$
$Q_k$ : Capacity of vehicle $k$
$M$ : A large constant (Big-M)

Decision Variables

$x_{ijk}$ : Binary variable, 1 if vehicle $k$ travels from node $i$ to node $j$ , 0 otherwise
$t_{ik}$ : Arrival time at node $i$ by vehicle $k$
$t_{start_k}$ : Departure time from depot for vehicle $k$
$t_{end_k}$ : Return time to depot for vehicle $k$
$y_{ik}$ : Binary variable, 1 if vehicle $k$ visits node $i$ , 0 otherwise
$u_{ik}$ : MTZ variable for subtour elimination

Objective Function

Minimize the total operation time of all vehicles:

\text{Minimize} \quad \sum_{k \in K} (t_{end_k} - t_{start_k})

Constraints

1. Each customer is visited exactly once

\sum_{k \in K} y_{ik} = 1 \quad \forall i \in V

2. Each vehicle departs from and returns to the depot exactly once

\sum_{j \in N, j \neq i} x_{0jk} = 1 \quad \forall k \in K

\sum_{i \in N, i \neq j} x_{i0k} = 1 \quad \forall k \in K

3. Flow conservation constraints for nodes (excluding depot)

\sum_{j \in N, j \neq i} x_{ijk} - \sum_{j \in N, j \neq i} x_{jik} = 0 \quad \forall i \in V, \forall k \in K

4. Link $x$ and $y$ variables

\sum_{j \in N, j \neq i} x_{ijk} = y_{ik} \quad \forall i \in V, \forall k \in K

5. Capacity constraints for each vehicle

\sum_{i \in V} d_i y_{ik} \leq Q_k \quad \forall k \in K

6. Time window constraints and time propagation constraints

Departure time from depot is within time window:

$a_0 \leq t_{start_k} \leq b_0 \quad \forall k \in K$
Time propagation constraints:
- From depot to first customer: $t_{jk} \geq t_{start_k} + c_{0j} - M(1 - x_{0jk}) \quad \forall j \in V, \forall k \in K$
- Between customers: $t_{jk} \geq t_{ik} + c_{ij} - M(1 - x_{ijk}) \quad \forall i, j \in V, i \neq j, \forall k \in K$
- Return to depot: $t_{end_k} \geq t_{ik} + c_{i0} - M(1 - x_{i0k}) \quad \forall i \in V, \forall k \in K$
Time window constraints for customer arrivals:
$a_i \leq t_{ik} \leq b_i \quad \forall i \in V, \forall k \in K$

7. Subtour elimination (MTZ) constraints

MTZ constraints for subtour elimination: $u_{ik} \geq 1 \quad \forall i \in V, \forall k \in K$ $u_{ik} \leq n \quad \forall i \in V, \forall k \in K$ $u_{ik} - u_{jk} + n x_{ijk} \leq n - 1 \quad \forall i, j \in V, i \neq j, \forall k \in K$

8. Ensure that if a customer is not visited by vehicle $k$ , its time variables are inactive

y_{ik} \leq \sum_{j \in N, j \neq i} x_{jik} \quad \forall i \in V, \forall k \in K

Lazy Constraint Generation for Subtour Elimination

Here is the mathematical formulation, I implemented in this code.

from pyscipopt import Model, quicksum, Conshdlr, SCIP_RESULT
import networkx as nx

def create_data_model():
    """Stores the data for the problem."""
    data = {}
    # Time matrix representing travel times between locations (consistent with time windows)
    data["time_matrix"] = [
        [0, 6, 9, 8, 7, 3, 6, 2, 3, 2, 6, 6, 4, 4, 5, 9, 7],
        [6, 0, 8, 3, 2, 6, 8, 4, 8, 8, 13, 7, 5, 8, 12, 10, 14],
        [9, 8, 0, 11, 10, 6, 3, 9, 5, 8, 4, 15, 14, 13, 9, 18, 9],
        [8, 3, 11, 0, 1, 7, 10, 6, 10, 10, 14, 6, 7, 9, 14, 6, 16],
        [7, 2, 10, 1, 0, 6, 9, 4, 8, 9, 13, 4, 6, 8, 12, 8, 14],
        [3, 6, 6, 7, 6, 0, 2, 3, 2, 2, 7, 9, 7, 7, 6, 12, 8],
        [6, 8, 3, 10, 9, 2, 0, 6, 2, 5, 4, 12, 10, 10, 6, 15, 5],
        [2, 4, 9, 6, 4, 3, 6, 0, 4, 4, 8, 5, 4, 3, 7, 8, 10],
        [3, 8, 5, 10, 8, 2, 2, 4, 0, 3, 4, 9, 8, 7, 3, 13, 6],
        [2, 8, 8, 10, 9, 2, 5, 4, 3, 0, 4, 6, 5, 4, 3, 9, 5],
        [6, 13, 4, 14, 13, 7, 4, 8, 4, 4, 0, 10, 9, 8, 4, 13, 4],
        [6, 7, 15, 6, 4, 9, 12, 5, 9, 6, 10, 0, 1, 3, 7, 3, 10],
        [4, 5, 14, 7, 6, 7, 10, 4, 8, 5, 9, 1, 0, 2, 6, 4, 8],
        [4, 8, 13, 9, 8, 7, 10, 3, 7, 4, 8, 3, 2, 0, 4, 5, 6],
        [5, 12, 9, 14, 12, 6, 6, 7, 3, 3, 4, 7, 6, 4, 0, 9, 2],
        [9, 10, 18, 6, 8, 12, 15, 8, 13, 9, 13, 3, 4, 5, 9, 0, 9],
        [7, 14, 9, 16, 14, 8, 5, 10, 6, 5, 4, 10, 8, 6, 2, 9, 0],
    ]
    data["time_windows"] = [
        (0, 5),     # 0 - 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
        (10, 16),   # 10
        (10, 15),   # 11
        (0, 5),     # 12
        (5, 10),    # 13
        (7, 8),     # 14
        (10, 15),   # 15
        (11, 15),   # 16
    ]
    data["demands"] = [0, 1, 1, 2, 4, 2, 4, 8, 8, 1, 2, 1, 2, 4, 4, 8, 8]
    data["vehicle_capacities"] = [15, 15, 15, 15]  # Capacities for each vehicle
    data["num_vehicles"] = 4
    data["depot"] = 0
    return data

def main():
    """Solve the VRPTW problem using PySCIPOpt with subtour elimination as lazy constraints."""
    # Retrieve data
    data = create_data_model()
    time_matrix = data["time_matrix"]
    time_windows = data["time_windows"]
    demands = data["demands"]
    vehicle_capacities = data["vehicle_capacities"]
    num_vehicles = data["num_vehicles"]
    depot = data["depot"]

    n = len(time_matrix)  # Number of nodes
    N = range(n)  # Set of nodes
    V = [i for i in N if i != depot]  # Customers excluding depot
    K = range(num_vehicles)  # Set of vehicles

    # Calculate M for Big-M constraints (maximum possible time)
    max_time = max(tw[1] for tw in time_windows) + max(
        time_matrix[i][j] for i in N for j in N
    )

    # Initialize the model
    model = Model("VRPTW with Lazy Constraints")

    # Variables
    x = {}  # Binary variables: x[i,j,k] == 1 if vehicle k travels from i to j
    t = {}  # Arrival time at node i by vehicle k
    t_start = {}  # Departure time from depot for vehicle k
    t_end = {}    # Return time to depot for vehicle k
    y = {}  # Binary variables: y[i,k] == 1 if vehicle k visits node i

    # Decision variables
    for k in K:
        capacity_k = vehicle_capacities[k]
        # Departure time from depot
        t_start[k] = model.addVar(lb=time_windows[depot][0], ub=time_windows[depot][1], vtype="I", name=f"t_start_{k}")
        # Return time to depot
        t_end[k] = model.addVar(lb=time_windows[depot][0], ub=None, vtype="I", name=f"t_end_{k}")
        for i in N:
            if i != depot:
                y[i, k] = model.addVar(vtype="B", name=f"y_{i}_{k}")
                t[i, k] = model.addVar(lb=time_windows[i][0], ub=time_windows[i][1], vtype="I", name=f"t_{i}_{k}")
            for j in N:
                if i != j:
                    x[i, j, k] = model.addVar(vtype="B", name=f"x_{i}_{j}_{k}")

    # Objective Function: Minimize total operation time of all vehicles
    model.setObjective(
        quicksum(t_end[k] - t_start[k] for k in K),
        "minimize",
    )

    # Constraints

    # 1. Each customer is visited exactly once
    for i in V:
        model.addCons(
            quicksum(y[i, k] for k in K) == 1,
            name=f"VisitOnce_{i}"
        )

    # 2. For each vehicle, ensure it departs from and returns to the depot correctly
    for k in K:
        model.addCons(
            quicksum(x[depot, j, k] for j in N if j != depot) == 1,
            name=f"DepartDepot_{k}"
        )
        model.addCons(
            quicksum(x[i, depot, k] for i in N if i != depot) == 1,
            name=f"ReturnDepot_{k}"
        )

    # 3. Flow conservation constraints for nodes (excluding depot)
    for k in K:
        for i in V:
            model.addCons(
                quicksum(x[i, j, k] for j in N if j != i) - quicksum(x[j, i, k] for j in N if j != i) == 0,
                name=f"FlowConservation_{i}_{k}"
            )

    # 4. Link x and y variables
    for k in K:
        for i in V:
            model.addCons(
                quicksum(x[i, j, k] for j in N if j != i) == y[i, k],
                name=f"LinkXY_{i}_{k}"
            )

    # 5. Capacity constraints for each vehicle
    for k in K:
        capacity_k = vehicle_capacities[k]
        model.addCons(
            quicksum(demands[i] * y[i, k] for i in V) <= capacity_k,
            name=f"VehicleCapacity_{k}"
        )

    # 6. Time window constraints and time propagation constraints
    for k in K:
        # Departure time from depot is within time window
        model.addCons(
            t_start[k] >= time_windows[depot][0],
            name=f"DepotTimeWindowStart_{k}"
        )
        model.addCons(
            t_start[k] <= time_windows[depot][1],
            name=f"DepotTimeWindowEnd_{k}"
        )

        # Time propagation constraints
        for i in N:
            if i != depot:
                # Arrival time at node i for vehicle k
                model.addCons(
                    t[i, k] >= time_windows[i][0],
                    name=f"TimeWindowStart_{i}_{k}"
                )
                model.addCons(
                    t[i, k] <= time_windows[i][1],
                    name=f"TimeWindowEnd_{i}_{k}"
                )
            for j in N:
                if i != j:
                    travel_time = time_matrix[i][j]
                    if i == depot:
                        # From depot to first customer
                        model.addCons(
                            t[j, k] >= t_start[k] + travel_time - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagationDepot_{i}_{j}_{k}"
                        )
                    elif j == depot:
                        # From last customer to depot
                        model.addCons(
                            t_end[k] >= t[i, k] + travel_time - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagationReturnDepot_{i}_{j}_{k}"
                        )
                    else:
                        # Between customers
                        model.addCons(
                            t[j, k] >= t[i, k] + travel_time - max_time * (1 - x[i, j, k]),
                            name=f"TimePropagation_{i}_{j}_{k}"
                        )

    # Include the subtour elimination constraint handler
    class SubtourEliminationConshdlr(Conshdlr):
        def __init__(self, x, n, N, V, K, depot):
            super().__init__()
            self.x = x
            self.n = n
            self.N = N
            self.V = V
            self.K = K
            self.depot = depot

        def conscheck(self, constraints, solution, checkintegrality, checklprows, printreason, completely):
            # For integer solutions, we ensure no subtours exist per vehicle
            for k in self.K:
                selected_edges = [(i, j) for i in self.N for j in self.N if i != j and solution[self.x[i, j, k]] > 0.5]
                G = nx.DiGraph()
                G.add_edges_from(selected_edges)
                # Find strongly connected components (subtours)
                subtours = list(nx.strongly_connected_components(G))
                for subtour in subtours:
                    if self.depot in subtour:
                        continue  # Skip component containing depot
                    if len(subtour) < 2:
                        continue
                    # Subtour detected, infeasible
                    if printreason:
                        print(f"Subtour detected for vehicle {k}: {subtour}")
                    return {"result": SCIP_RESULT.INFEASIBLE}

            return {"result": SCIP_RESULT.FEASIBLE}

        def consenfolp(self, constraints, nusefulconss, solinfeasible):
            # Enforce LP solution
            violated = False
            for k in self.K:
                selected_edges = [(i, j) for i in self.N for j in self.N if i != j and self.model.getVal(self.x[i, j, k]) > 0.5]
                G = nx.DiGraph()
                G.add_edges_from(selected_edges)
                # Find strongly connected components (subtours)
                subtours = list(nx.strongly_connected_components(G))
                for subtour in subtours:
                    if self.depot in subtour:
                        continue  # Skip component containing depot
                    if len(subtour) < 2:
                        continue
                    # Subtour detected, add subtour elimination constraint
                    subtour_nodes = list(subtour)
                    expr = quicksum(self.x[i, j, k] for i in subtour_nodes for j in subtour_nodes if i != j)
                    self.model.addCons(expr <= len(subtour_nodes) - 1)
                    violated = True
            if violated:
                return {"result": SCIP_RESULT.CONSADDED}
            else:
                return {"result": SCIP_RESULT.FEASIBLE}

        def consenfeas(self, constraints, solution, objinfeasible):
            # Enforce constraints on feasible integer solutions
            violated = False
            for k in self.K:
                selected_edges = [(i, j) for i in self.N for j in self.N if i != j and self.model.getSolVal(solution, self.x[i, j, k]) > 0.5]
                G = nx.DiGraph()
                G.add_edges_from(selected_edges)
                # Find strongly connected components (subtours)
                subtours = list(nx.strongly_connected_components(G))
                for subtour in subtours:
                    if self.depot in subtour:
                        continue  # Skip component containing depot
                    if len(subtour) < 2:
                        continue
                    # Subtour detected, add subtour elimination constraint
                    subtour_nodes = list(subtour)
                    expr = quicksum(self.x[i, j, k] for i in subtour_nodes for j in subtour_nodes if i != j)
                    self.model.addCons(expr <= len(subtour_nodes) - 1)
                    violated = True
            if violated:
                return {"result": SCIP_RESULT.CONSADDED}
            else:
                return {"result": SCIP_RESULT.FEASIBLE}

    # Include the constraint handler
    conshdlr = SubtourEliminationConshdlr(x, n, N, V, K, depot)
    model.includeConshdlr(conshdlr, "SubtourEliminationConshdlr", "Eliminates subtours in VRPTW")

    # Solve the model
    model.optimize()

    # Retrieve and print the solution
    status = model.getStatus()
    if status == "optimal" or status == "bestsollimit":
        print_solution(data, x, t_start, t_end, t, y, K, N, depot, model)
    else:
        print("No optimal solution found.")
        print(f"Solver status: {status}")

def print_solution(data, x, t_start, t_end, t, y, K, N, depot, model):
    """Prints solution in the specified format."""
    solution = model.getBestSol()
    total_time = 0

    print(f"Objective: {model.getSolObjVal(solution):.1f}")
    for k in K:
        vehicle_id = k
        # Build the route for vehicle k
        edges = [(i, j) for i in N for j in N if i != j and model.getSolVal(solution, x.get((i, j, k), 0)) > 0.5]
        if not edges:
            continue  # Skip if vehicle k is not used
        G = nx.DiGraph()
        G.add_edges_from(edges)

        if depot not in G.nodes:
            continue  # Skip if vehicle k is not used

        plan_output = f"Route for vehicle {vehicle_id}:\n"

        route = [depot]
        time_log = [f"Time({model.getSolVal(solution, t_start[k]):.1f})"]
        next_node = depot

        while True:
            successors = list(G.successors(next_node))
            if not successors:
                break
            next_node = successors[0]
            if next_node == depot:
                break
            route.append(next_node)
            arrival_time = model.getSolVal(solution, t[next_node, k])
            time_log.append(f"Time({arrival_time:.1f})")

        # Append depot at the end
        route.append(depot)
        time_log.append(f"Time({model.getSolVal(solution, t_end[k]):.1f})")

        # Print the route and times
        for idx, node in enumerate(route):
            plan_output += f"{node} {time_log[idx]} -> "
        plan_output = plan_output.rstrip(" -> ") + "\n"
        route_duration = model.getSolVal(solution, t_end[k]) - model.getSolVal(solution, t_start[k])
        plan_output += f"Time of the route: {route_duration:.1f}min\n"
        print(plan_output)
        total_time += route_duration

    print(f"Total time of all routes: {total_time:.1f}min")

if __name__ == "__main__":
    main()

Sets and Indices

$N$ : Set of all nodes (including depot)
$V$ : Set of customer nodes (excluding depot)
$K$ : Set of vehicles
$i, j \in N$ : Indices for nodes
$k \in K$ : Index for vehicles

Parameters

$c_{ij}$ : Travel time from node $i$ to node $j$
$[a_i, b_i]$ : Time window for node $i$
$d_i$ : Demand at node $i$
$Q_k$ : Capacity of vehicle $k$
$M$ : A large constant (Big-M)

Decision Variables

$x_{ijk}$ : Binary variable, 1 if vehicle $k$ travels from node $i$ to node $j$ , 0 otherwise
$t_{ik}$ : Arrival time at node $i$ by vehicle $k$
$t_{start_k}$ : Departure time from depot for vehicle $k$
$t_{end_k}$ : Return time to depot for vehicle $k$
$y_{ik}$ : Binary variable, 1 if vehicle $k$ visits node $i$ , 0 otherwise

Objective Function

Minimize the total operation time of all vehicles:

\text{Minimize} \quad \sum_{k \in K} (t_{end_k} - t_{start_k})

Constraints

1. Each customer is visited exactly once

\sum_{k \in K} y_{ik} = 1 \quad \forall i \in V

2. Each vehicle departs from and returns to the depot exactly once

\sum_{j \in N, j \neq i} x_{0jk} = 1 \quad \forall k \in K

\sum_{i \in N, i \neq j} x_{i0k} = 1 \quad \forall k \in K

3. Flow conservation constraints for nodes (excluding depot)

\sum_{j \in N, j \neq i} x_{ijk} - \sum_{j \in N, j \neq i} x_{jik} = 0 \quad \forall i \in V, \forall k \in K

4. Link $x$ and $y$ variables

\sum_{j \in N, j \neq i} x_{ijk} = y_{ik} \quad \forall i \in V, \forall k \in K

5. Capacity constraints for each vehicle

\sum_{i \in V} d_i y_{ik} \leq Q_k \quad \forall k \in K

6. Time window constraints and time propagation constraints

Departure time from depot is within time window:

$a_0 \leq t_{start_k} \leq b_0 \quad \forall k \in K$
Time propagation constraints:
- From depot to first customer: $t_{jk} \geq t_{start_k} + c_{0j} - M(1 - x_{0jk}) \quad \forall j \in V, \forall k \in K$
- Between customers: $t_{jk} \geq t_{ik} + c_{ij} - M(1 - x_{ijk}) \quad \forall i, j \in V, i \neq j, \forall k \in K$
- Return to depot: $t_{end_k} \geq t_{ik} + c_{i0} - M(1 - x_{i0k}) \quad \forall i \in V, \forall k \in K$
Time window constraints for customer arrivals:
$a_i \leq t_{ik} \leq b_i \quad \forall i \in V, \forall k \in K$

7. Subtour elimination constraints (Lazy Constraint Generation)

Subtour elimination constraints are dynamically added during the optimization process. For any detected subtour $S$ : $\sum_{i \in S} \sum_{j \in S, j \neq i} x_{ijk} \leq |S| - 1 \quad \forall k \in K$

Dynamic Addition of Constraints

In this model, subtour elimination constraints are not added upfront but are dynamically generated during the optimization process. This approach is known as Lazy Constraint Generation. During the optimization, the solver periodically checks the current solution for any subtours. If a subtour is detected, a new constraint is added to the model to eliminate that subtour. This process continues until no more subtours are found, ensuring that the final solution is feasible and optimal.

Performance Comparison

Below is a table comparing the performance of different models and tools used to solve the Capacitated Vehicle Routing Problem with Time Windows (CVRPTW). The comparison includes execution time, lines of code, and the total time of all routes.

Model/Tool	Execution Time in secs	Lines of Code	Total Time of All Routes in mins (Objective Function Value)
OR-Tools	0.02	170	82
Node-based Model (SCIP)	0.54	283	81
Flow-based Model (SCIP)	21	294	81
MTZ Model (SCIP)	11	288	81
Lazy Constraint Generation for Subtour Elimination (SCIP)	4	339	81

Description

OR-Tools: The fastest among all, with an execution time of 0.02 seconds and 170 lines of code. However, the total time of all routes is slightly higher at 82 minutes.
Node-based Model (SCIP): Provides a slightly better solution with a total route time of 81 minutes and executes in 0.54 seconds. The implementation requires 283 lines of code.
Flow-based Model (SCIP): Takes significantly longer to execute at 21 seconds but achieves the same total route time of 81 minutes. The implementation consists of 294 lines of code.
MTZ Model (SCIP): Executes in 11 seconds with a total route time of 81 minutes and requires 288 lines of code.
Lazy Constraint Generation for Subtour Elimination (SCIP): Balances execution time and code complexity, taking 4 seconds to execute with 339 lines of code, and achieves a total route time of 81 minutes.

This comparison highlights the trade-offs between execution time, code complexity, and solution quality. OR-Tools is the fastest but slightly less optimal in terms of total route time. SCIP models, particularly the Node-based model, provide better solutions but at the cost of increased execution time and code complexity.

Personal Experience and Insights

Ease of Use and Code Complexity

OR-Tools

Fewer Lines of Code: One of the standout features of OR-Tools is its simplicity and ease of use. The API is designed to be user-friendly, which allows for rapid development and fewer lines of code. This makes OR-Tools an excellent choice for proof-of-concept (PoC) projects or quick benchmarks.
Good for PoC or Benchmark: Due to its simplicity, OR-Tools is ideal for quickly setting up and testing optimization problems. It allows you to focus on the problem at hand without getting bogged down by the intricacies of the underlying algorithms.

SCIP

More Control and Customization: SCIP offers a high degree of control and customization, allowing you to fine-tune the optimization process to meet specific requirements. This is particularly useful for complex problems where standard solvers may not perform well.
More Complex: The trade-off for this level of control is increased complexity. SCIP requires a deeper understanding of optimization techniques and more lines of code to implement. This can be a barrier for beginners but is invaluable for advanced users who need to customize their solutions.

Reliability and Support

OR-Tools

Risk of Google Discontinuing Projects: One of the potential drawbacks of using OR-Tools is the risk associated with Google's history of discontinuing projects. While OR-Tools is currently well-supported and actively maintained, there is always a possibility that Google may decide to discontinue the project in the future. This uncertainty can be a concern for users who rely on OR-Tools for critical applications.
Bad Documentation: Another challenge with OR-Tools is the quality of its documentation. While the official documentation provides a good starting point, it can be lacking in detail and clarity, making it difficult for users to fully understand and utilize all the features and capabilities of the tool. This can be particularly frustrating for new users who are trying to learn how to use OR-Tools effectively.
Bad Choice for Deployment: Given the potential risks and challenges associated with OR-Tools, it may not be the best choice for deployment in production environments. The uncertainty around long-term support and the limitations of the documentation can make it difficult to ensure the reliability and maintainability of applications built with OR-Tools.

SCIP

Good Choice for Deployment: In contrast, SCIP is a more reliable choice for deployment in production environments. As an open-source project with a strong academic and industrial user base, SCIP is less likely to be discontinued or lose support. Additionally, the comprehensive documentation and active community make it easier for users to get help and find solutions to any issues they may encounter.
Better Documentation: SCIP provides extensive and detailed documentation, including examples and tutorials, which can help users understand and effectively utilize the tool. This makes it easier to implement and maintain optimization models, particularly for complex and specialized problems.
Long-Term Support: The strong academic and industrial backing of SCIP ensures that it will continue to be supported and developed for the foreseeable future. This makes it a more reliable choice for long-term projects and critical applications where stability and support are essential.

Conclusion

In summary, while OR-Tools offers simplicity and ease of use, making it a good choice for quick proofs of concept and benchmarking, it may not be the best option for deployment due to potential risks and documentation challenges. On the other hand, SCIP provides greater control and customization, along with better documentation and long-term support, making it a more reliable choice for deployment in production environments.

Introduction

Problem Statement

Tools Overview

Google OR-Tools

Brief Introduction to OR-Tools

Key Features and Advantages

Use Cases and Applications

SCIP

Brief Introduction to SCIP

Key Features and Advantages

Use Cases and Applications

Implementation Details

Google OR-Tools Implementation

Routing Index Manager

Creating the Routing Model

Registering the Time Callback

Setting the Objective Function

Registering the Demand Callback

Adding the Capacity Constraint

Adding the Time Windows Constraint

Adding Time Window Constraints for Each Location

Adding Time Window Constraints for Each Vehicle Start Node

Instantiating Route Start and End Times

Optimizing

SCIP Implementations

Node-Based Model

Sets and Indices

Parameters

Decision Variables

Objective Function

Constraints

1. Each customer is visited exactly once

2. Each vehicle departs from and returns to the depot exactly once

3. Link xxx and yyy variables

4. Flow conservation constraints for nodes (excluding depot)

5. Capacity constraints for each vehicle

6. Time window constraints and time propagation constraints

Departure time from depot is within time window

Time propagation constraints

Flow-based Model

Sets and Indices

Parameters

Decision Variables

Objective Function

Constraints

1. Each customer is visited exactly once

2. Each vehicle departs from and returns to the depot exactly once

3. Link xxx and yyy variables

4. Flow conservation constraints (using flow variables)

5. Time window constraints and time propagation constraints

6. Time windows for customer arrivals

MTZ Model

Sets and Indices

Parameters

Decision Variables

Objective Function

Constraints

1. Each customer is visited exactly once

2. Each vehicle departs from and returns to the depot exactly once

3. Flow conservation constraints for nodes (excluding depot)

4. Link xxx and yyy variables

5. Capacity constraints for each vehicle

6. Time window constraints and time propagation constraints

7. Subtour elimination (MTZ) constraints

8. Ensure that if a customer is not visited by vehicle kkk , its time variables are inactive

Lazy Constraint Generation for Subtour Elimination

Sets and Indices

Parameters

Decision Variables

Objective Function

Constraints

1. Each customer is visited exactly once

2. Each vehicle departs from and returns to the depot exactly once

3. Flow conservation constraints for nodes (excluding depot)

4. Link xxx and yyy variables

5. Capacity constraints for each vehicle

6. Time window constraints and time propagation constraints

7. Subtour elimination constraints (Lazy Constraint Generation)

Dynamic Addition of Constraints

Performance Comparison

3. Link $x$ and $y$ variables

3. Link $x$ and $y$ variables

4. Link $x$ and $y$ variables

8. Ensure that if a customer is not visited by vehicle $k$ , its time variables are inactive

4. Link $x$ and $y$ variables