I have a question with multiple locations, time windows for each location, multiple vehicles, and weight capacity for each vehicle. I need to get a vehicle routing result for each vehicle considering the capacity and time windows of locations.

I started with the or-tool example cvrptw.py. and try to add distance matrix as input using the or-tool example vrp_capacity.py. I added the changed code below.

from functools import partial
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['distance_matrix'] = [
    [
        0, 548, 776, 696, 582, 274, 502, 194, 308, 194, 536, 502, 388, 354,
        468, 776, 662
    ],
    [
        548, 0, 684, 308, 194, 502, 730, 354, 696, 742, 1084, 594, 480, 674,
        1016, 868, 1210
    ],
    [
        776, 684, 0, 992, 878, 502, 274, 810, 468, 742, 400, 1278, 1164,
        1130, 788, 1552, 754
    ],
    [
        696, 308, 992, 0, 114, 650, 878, 502, 844, 890, 1232, 514, 628, 822,
        1164, 560, 1358
    ],
    [
        582, 194, 878, 114, 0, 536, 764, 388, 730, 776, 1118, 400, 514, 708,
        1050, 674, 1244
    ],
    [
        274, 502, 502, 650, 536, 0, 228, 308, 194, 240, 582, 776, 662, 628,
        514, 1050, 708
    ],
    [
        502, 730, 274, 878, 764, 228, 0, 536, 194, 468, 354, 1004, 890, 856,
        514, 1278, 480
    ],
    [
        194, 354, 810, 502, 388, 308, 536, 0, 342, 388, 730, 468, 354, 320,
        662, 742, 856
    ],
    [
        308, 696, 468, 844, 730, 194, 194, 342, 0, 274, 388, 810, 696, 662,
        320, 1084, 514
    ],
    [
        194, 742, 742, 890, 776, 240, 468, 388, 274, 0, 342, 536, 422, 388,
        274, 810, 468
    ],
    [
        536, 1084, 400, 1232, 1118, 582, 354, 730, 388, 342, 0, 878, 764,
        730, 388, 1152, 354
    ],
    [
        502, 594, 1278, 514, 400, 776, 1004, 468, 810, 536, 878, 0, 114,
        308, 650, 274, 844
    ],
    [
        388, 480, 1164, 628, 514, 662, 890, 354, 696, 422, 764, 114, 0, 194,
        536, 388, 730
    ],
    [
        354, 674, 1130, 822, 708, 628, 856, 320, 662, 388, 730, 308, 194, 0,
        342, 422, 536
    ],
    [
        468, 1016, 788, 1164, 1050, 514, 514, 662, 320, 274, 388, 650, 536,
        342, 0, 764, 194
    ],
    [
        776, 868, 1552, 560, 674, 1050, 1278, 742, 1084, 810, 1152, 274,
        388, 422, 764, 0, 798
    ],
    [
        662, 1210, 754, 1358, 1244, 708, 480, 856, 514, 468, 354, 844, 730,
        536, 194, 798, 0
    ],
]
data['num_locations'] = len(data['distance_matrix'])
data['time_windows'] = \
      [(0, 0),
       (75, 85), (75, 85), 
       (60, 70), (45, 55), 
       (0, 8), (50, 60),
       (0, 10), (10, 20), 
       (0, 10), (75, 85), 
       (85, 95), (5, 15), 
       (15, 25), (10, 20), 
       (45, 55), (30, 40)] 
data['demands'] = \
      [0, 
       1, 1, 
       2, 4, 
       2, 4, 
       8, 8, 
       1, 2, 
       1, 2, 
       4, 4, 
       8, 8] 
data['time_per_demand_unit'] = 5  
data['num_vehicles'] = 4
data['vehicle_capacity'] = 15
data['vehicle_speed'] = 83  
data['depot'] = 0
return data




def manhattan_distance(position_1, position_2):

return (
    abs(position_1[0] - position_2[0]) + abs(position_1[1] - position_2[1]))





def create_demand_evaluator(data):

_demands = data['demands']

def demand_evaluator(manager, node):
   
    return _demands[manager.IndexToNode(node)]

return demand_evaluator


def add_capacity_constraints(routing, data, demand_evaluator_index):

capacity = 'Capacity'
routing.AddDimension(
    demand_evaluator_index,
    0,  # null capacity slack
    data['vehicle_capacity'],
    True,  # start cumul to zero
    capacity)


def create_time_evaluator(data):

def service_time(data, node):
   
    return data['demands'][node] * data['time_per_demand_unit']

def travel_time(data, from_node, to_node):
    """Gets the travel times between two locations."""
    if from_node == to_node:
        travel_time = 0
    else:
        travel_time = manhattan_distance(data['locations'][from_node], data[
            'locations'][to_node]) / data['vehicle_speed']
    return travel_time

_total_time = {}
# precompute total time to have time callback in O(1)
for from_node in range(data['num_locations']):
    _total_time[from_node] = {}
    for to_node in range(data['num_locations']):
        if from_node == to_node:
            _total_time[from_node][to_node] = 0
        else:
            _total_time[from_node][to_node] = int(
                service_time(data, from_node) + travel_time(
                    data, from_node, to_node))

def time_evaluator(manager, from_node, to_node):
    """Returns the total time between the two nodes"""
    return _total_time[manager.IndexToNode(from_node)][manager.IndexToNode(
        to_node)]

return time_evaluator


def add_time_window_constraints(routing, manager, data, time_evaluator_index):
"""Add Global Span constraint"""
time = 'Time'
horizon = 120
routing.AddDimension(
    time_evaluator_index,
    horizon,  # allow waiting time
    horizon,  # maximum time per vehicle
    False,  # don't force start cumul to zero since we are giving TW to start nodes
    time)
time_dimension = routing.GetDimensionOrDie(time)
# Add time window constraints for each location except depot
# and 'copy' the slack var in the solution object (aka Assignment) to print it
for location_idx, time_window in enumerate(data['time_windows']):
    if location_idx == 0:
        continue
    index = manager.NodeToIndex(location_idx)
    time_dimension.CumulVar(index).SetRange(time_window[0], time_window[1])
    routing.AddToAssignment(time_dimension.SlackVar(index))
# Add time window constraints for each vehicle start node
# and 'copy' the slack var in the solution object (aka Assignment) to print it
for vehicle_id in range(data['num_vehicles']):
    index = routing.Start(vehicle_id)
    time_dimension.CumulVar(index).SetRange(data['time_windows'][0][0],
                                            data['time_windows'][0][1])
    routing.AddToAssignment(time_dimension.SlackVar(index))
    # Warning: Slack var is not defined for vehicle's end node
    #routing.AddToAssignment(time_dimension.SlackVar(self.routing.End(vehicle_id)))



def print_solution(manager, routing, assignment):  # pylint:disable=too-many-locals
"""Prints assignment on console"""
print(f'Objective: {assignment.ObjectiveValue()}')
time_dimension = routing.GetDimensionOrDie('Time')
capacity_dimension = routing.GetDimensionOrDie('Capacity')
total_distance = 0
total_load = 0
total_time = 0
for vehicle_id in range(manager.GetNumberOfVehicles()):
    index = routing.Start(vehicle_id)
    plan_output = 'Route for vehicle {}:\n'.format(vehicle_id)
    distance = 0
    while not routing.IsEnd(index):
        load_var = capacity_dimension.CumulVar(index)
        time_var = time_dimension.CumulVar(index)
        slack_var = time_dimension.SlackVar(index)
        plan_output += ' {0} Load({1}) Time({2},{3}) Slack({4},{5}) ->'.format(
            manager.IndexToNode(index),
            assignment.Value(load_var),
            assignment.Min(time_var),
            assignment.Max(time_var),
            assignment.Min(slack_var), assignment.Max(slack_var))
        previous_index = index
        index = assignment.Value(routing.NextVar(index))
        distance += routing.GetArcCostForVehicle(previous_index, index,
                                                 vehicle_id)
    load_var = capacity_dimension.CumulVar(index)
    time_var = time_dimension.CumulVar(index)
    slack_var = time_dimension.SlackVar(index)
    plan_output += ' {0} Load({1}) Time({2},{3})\n'.format(
        manager.IndexToNode(index),
        assignment.Value(load_var),
        assignment.Min(time_var), assignment.Max(time_var))
    plan_output += 'Distance of the route: {0}m\n'.format(distance)
    plan_output += 'Load of the route: {}\n'.format(
        assignment.Value(load_var))
    plan_output += 'Time of the route: {}\n'.format(
        assignment.Value(time_var))
    print(plan_output)
    total_distance += distance
    total_load += assignment.Value(load_var)
    total_time += assignment.Value(time_var)
print('Total Distance of all routes: {0}m'.format(total_distance))
print('Total Load of all routes: {}'.format(total_load))
print('Total Time of all routes: {0}min'.format(total_time))
# [END solution_printer]


def main():
"""Solve the Capacitated VRP with time windows."""
# Instantiate the data problem.
# [START data]
data = create_data_model()
# [END data]

# Create the routing index manager.
# [START index_manager]
manager = pywrapcp.RoutingIndexManager(data['num_locations'],
                                       data['num_vehicles'], data['depot'])
# [END index_manager]

# Create Routing Model.
# [START routing_model]
routing = pywrapcp.RoutingModel(manager)
# [END routing_model]

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)
# [END transit_callback]

# Define cost of each arc.
# [START arc_cost]
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)


# Add Capacity constraint.
# [START capacity_constraint]
demand_evaluator_index = routing.RegisterUnaryTransitCallback(
    partial(create_demand_evaluator(data), manager))
add_capacity_constraints(routing, data, demand_evaluator_index)
# [END capacity_constraint]

# Add Time Window constraint.
# [START time_constraint]
time_evaluator_index = routing.RegisterTransitCallback(
    partial(create_time_evaluator(data), manager))
add_time_window_constraints(routing, manager, data, time_evaluator_index)
# [END time_constraint]

# Setting first solution heuristic (cheapest addition).
# [START parameters]
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(2)
search_parameters.log_search = True
# [END parameters]

# Solve the problem.
# [START solve]
solution = routing.SolveWithParameters(search_parameters)
# [END solve]

# Print solution on console.
# [START print_solution]
if solution:
    print_solution(manager, routing, solution)
else:
    print('No solution found!')
# [END print_solution]


if __name__ == '__main__':
main()

I changed,

1. Replace _locations array with distance matrix array
2. Remove distance_evaluator_index and add distance callback.

Now I need to modify create_time_evaluator according to the distance matrix and apply it to time_evaluator_index. How can I achieve this and get better results?

I need to achieve Capacity constraints and time window constraints using distance matrix as a single solution with minimum cost. Is this the better way for my problem?