CVRP with delivery time

[tcokyasar44]

Cross-posted

I am modeling a CVRP with a single vehicle type. The objective is to minimize the travel time. Total time (travel + delivery/service at nodes) is capped by 36,000 for each vehicle. Each node demands 3 units of deliveries, and the vehicle capacity is 120 units for each vehicle.

I programmed it in Python following the CVRP examples. The code is below. As a newbie, I have the below listed questions:

1. Time_matrix data seems to be converted to integers though I use floats. Is it possible to change this?
2. How can I add a node-dependent delivery/service time that does not count on to the objective function? In other words, I would like the service time to be counted towards the total time constraint but not in the objective.
3. When I use a time limit of 120 (as seen below), I got a solution that violates the vehicle capacity constraint by 1 unit, i.e., one of the vehicles delivered 121 packages. Is it possible to get a feasible solution after the time limit is reached?
4. Is it possible to receive a quick, (meta-) heuristic, feasible solution?
5. Is there a way to consider the number of vehicles a variable in my problem setting?
6. Is there an option to see the optimality gap (or at least gap from the lower bound)?
7. Is there a solver selection option, like Gurobi, CPLEX, etc.?
8. Is there a way to see the logs while solving?
9. Up to how many nodes should I expect a good solution in this setting? I understand it depends on the data, and so on; but, I am just chasing an estimate.

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 = {}
    TT = np.array(df.values)
    data['time_matrix'] = TT
    data['num_vehicles'] = 15
    data['depot'] = 0
    data['demands'] = [1]*len(TT)
    data['vehicle_capacities'] = [120]*data['num_vehicles']
    return data

def print_solution(data, manager, routing, solution):
    """Prints solution on console."""
    print(f'Objective: {solution.ObjectiveValue()}')
    max_route_time = 0; tour = {i:[] for i in range(data['num_vehicles'])}
    for vehicle_id in range(data['num_vehicles']):
        index = routing.Start(vehicle_id)
        plan_output = 'Route for vehicle {}:\n'.format(vehicle_id)
        route_time = 0
        while not routing.IsEnd(index):
            plan_output += ' {} -> '.format(manager.IndexToNode(index))
            tour[vehicle_id].append(manager.IndexToNode(index))
            previous_index = index
            index = solution.Value(routing.NextVar(index))
            route_time += routing.GetArcCostForVehicle(
                previous_index, index, vehicle_id)
        plan_output += '{}\n'.format(manager.IndexToNode(index))
        tour[vehicle_id].append(manager.IndexToNode(index))
        plan_output += 'Travel time of the route: {} sec\n'.format(route_time)
        print(plan_output)
        max_route_time = max(route_time, max_route_time)
    print('Maximum of the route time: {} sec'.format(max_route_time))
    return(tour)

"""Entry point of the program."""
# 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 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 Distance constraint.
dimension_name = 'Time'
routing.AddDimension(
    transit_callback_index,
    0,  # no slack
    36000,  # vehicle maximum travel time
    True,  # start cumul to zero
    dimension_name)
distance_dimension = routing.GetDimensionOrDie(dimension_name)

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

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

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

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

# Print solution on console.
if solution:
    tour = print_solution(data, manager, routing, solution)
else:
    print('No solution found !')

[Mizux]

1. Solver only support int64 so you must convert your float to int64 to avoid unintended behaviour/random crash
   note: you can multiply all your values by a coefficient e.g. *100 to simulate a 2 decimal precision before rounding etc...

2. By default, The objective will be only composed of the sum of all arcs cost used.
   So If I were you, I would create an evaluator without the service time to be use when setting the arc cost AND create a Time dimension with an other evaluator function for the time constraint with the service time.
   note: by default span coefficient of dimension is 0 i.e. arc of any AddDimension() won't be part of the objective function i.e. they will only be use to model constraint like maximum travel distance, Time Windows etc...

3. Strange, I'll try to investigate it, this is a first time I heard about this "bug"
   if you can provide a working example (i.e. with data) so we can reproduce it, it will be more than welcome ;)

4. No, but you can first run the solver with the GUIDED_LOCAL_SEARCH not enabled so solver will return once the first solution is found.
   note: Then you can use this solution as initial solution when performing a second run with this time the GUIDED_LOCAL_SEARCH enabled...

5. You can add a VehicleFixedCost so you set a big number of vehicle and see how solver try to minimize the objective by using the minimum vehicle as possible i.e. avoiding to pay too much fixed cost...

6. Not I'm aware of

7. Yes and No:
   No, the routing solver use a custom/ad_hoc constraint solver (that's why the routing source code is located in the constraint_solver directory ;))
   Yes, there is some on going work to use the CP-SAT solver to "finalize" the routing solution (especially for breaks and scheduling constraint management)
   see:

   or-tools/ortools/constraint_solver/routing_parameters.proto

   Lines 417 to 432 in b37d9c7

   	// If true, use the CP solver to find a solution. Either local or depth-first
   	// search will be used depending on the value of use_depth_first_search. Will
   	// be run before the CP-SAT solver (cf. use_cp_sat).
   	OptionalBoolean use_cp = 28;
   	// If true, use the CP-SAT solver to find a solution. If use_cp is also true,
   	// the CP-SAT solver will be run after the CP solver if there is time
   	// remaining and will use the CP solution as a hint for the CP-SAT search.
   	// As of 5/2019, only TSP models can be solved.
   	OptionalBoolean use_cp_sat = 27;
   	// Underlying solver to use in dimension scheduling, respectively for
   	// continuous and mixed models.
   	enum SchedulingSolver {
   	UNSET = 0;
   	GLOP = 1;
   	CP_SAT = 2;
   	}

8. yes

• In the routing level (very low level)

routing_parameters = pywrapcp.DefaultRoutingModelParameters()
routing_parameters.solver_parameters.trace_propagation = True
routing_parameters.solver_parameters.trace_search = True
routing = pywrapcp.RoutingModel(manager, routing_parameters)

• In the search parameter, high level, basically stats each time the GLS find a new solution

search_parameters.log_search = True
solution = routing.SolveWithParameters(search_parameters)

9. Don't get it, sorry
   by the way if all your nodes are not part of any Disjunction then it is mandatory for the solver to visit them all or return solution not found...
   If you want to allow solver to be able to drop nodes modulo paying a (high) penalty cost please take a look at example:
   https://developers.google.com/optimization/routing/penalties

[tcokyasar44]

1. 2 decimal is a satisfactory solution! Wouldn't it further increase the computational burden though?
2. If not too hard and time-consuming, a code snippet would be more helpful for this newbie.
3. It was my bad, I was counting the depot in. So, it was feasible.
4. Again, a code snippet would be appreciated.
5. I see the point, but I don't want to consider the cost of vehicles. I was actually wondering if making it a variable would lower the complexity. Wouldn't setting it to a large number increase complexity?
6. Then, this tool does not totally replace other optimization tools for the ones interested in the gap in routing problems.
7. Fair enough!
8. This is what I am talking about!
9. I am trying to solve instances with more than 1,500 nodes. Should I expect close to optimal results or God knows?

[tcokyasar44]

Let data['service_time'] denote the list of node-dependent service times. I think the below is the resolution to the node-dependent service time addition. Can you please double-check?

# Create and register a transit callback.
def time_callback(from_index, to_index):
    """Returns the 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)

# Create and register a transit callback.
def time_callback2(from_index, to_index):
    """Returns the 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)
    if from_node != 0:
        return data['time_matrix'][from_node][to_node] + data['service_time'][from_node]
    else:
        return data['time_matrix'][from_node][to_node]

transit_callback_index2 = routing.RegisterTransitCallback(time_callback2)

# Add Time constraint.
dimension_name = 'Time'
routing.AddDimension(
    transit_callback_index2,
    0,  # no slack
    6000,  # vehicle maximum travel time
    True,  # start cumul to zero
    dimension_name)
time_dimension = routing.GetDimensionOrDie(dimension_name)