Or-tools: Solution slack variables not working as expected

[ivarbrek]

OR tools v7.7 in Python, using the Routing Solver
OS: macOS Catalina 10.15.5

In the attached script I have modified cvrptw_break.py: mutate_locations make small random changes to the locations list, and number of vehicles are set to 6 to ensure feasibility. check_solution_output iterates every route to check that the equation slack(i) = cumul(j) - cumul(i) - transit(i, j) from the documentation is satisfied (with 1 time unit slack allowed):

def check_solution_output(data, manager, routing, assignment):
    time_dimension = routing.GetDimensionOrDie("Time")

    for vehicle_id in range(data["num_vehicles"]):
        i = routing.Start(vehicle_id)
        j = assignment.Value(routing.NextVar(i))

        while not routing.IsEnd(j):
            time_i = assignment.Value(time_dimension.CumulVar(i))
            time_j = assignment.Value(time_dimension.CumulVar(j))
            slack_i = assignment.Value(time_dimension.SlackVar(i))

            time_eval = create_time_evaluator(data)
            time_ij = time_eval(manager, i, j)

            assert time_j + 1 >= time_i + time_ij + slack_i >= time_j - 1, \
                "Time equation not in balance between nodes " + str(
                    manager.IndexToNode(i)) + "–" + str(manager.IndexToNode(j)) + ":  " + str(
                    time_j) + "!=" + str(time_i) + "+" + str(
                    time_ij) + "+" + str(slack_i) + "   indexes: " + str(i) + "-" + str(j)
            i = j
            j = assignment.Value(routing.NextVar(j))

I expected the provided solution to satisfy this condition, however, most times it does not. My output is typically: AssertionError: Time equation not in balance between nodes 5–4: 45!=3+13+24

From other cases this error is not present when slack in the time dimension is prohibited.

Is this a bug, or are there errors in my code?

Full code:

"""Capacitated Vehicle Routing Problem with Time Windows (CVRPTW).
   This is a sample using the routing library python wrapper to solve a CVRPTW
   problem.
   A description of the problem can be found here:
   http://en.wikipedia.org/wiki/Vehicle_routing_problem.
   Distances are in meters and time in minutes.
"""

from __future__ import print_function

from functools import partial

from ortools.constraint_solver import pywrapcp
from ortools.constraint_solver import routing_enums_pb2
from random import random, randint

def mutate_locations(locations):
    print(locations)
    for i, tup in enumerate(locations):
        new = list(tup[:])
        if random() < 0.2:
            new[0] += randint(-2, 2)
        if random() < 0.2:
            new[1] += randint(-2, 2)
        locations[i] = tuple(new)
    print(locations)
    return locations

def check_solution_output(data, manager, routing, assignment):
    time_dimension = routing.GetDimensionOrDie("Time")

    for vehicle_id in range(data["num_vehicles"]):
        i = routing.Start(vehicle_id)
        j = assignment.Value(routing.NextVar(i))

        while not routing.IsEnd(j):
            time_i = assignment.Value(time_dimension.CumulVar(i))
            time_j = assignment.Value(time_dimension.CumulVar(j))
            slack_i = assignment.Value(time_dimension.SlackVar(i))

            time_eval = create_time_evaluator(data)
            time_ij = time_eval(manager, i, j)

            assert time_j + 1 >= time_i + time_ij + slack_i >= time_j - 1, \
                "Time equation not in balance between nodes " + str(
                    manager.IndexToNode(i)) + "–" + str(manager.IndexToNode(j)) + ":  " + str(
                    time_j) + "!=" + str(time_i) + "+" + str(
                    time_ij) + "+" + str(slack_i) + "   indexes: " + str(i) + "-" + str(j)
            i = j
            j = assignment.Value(routing.NextVar(j))

###########################
# Problem Data Definition #
###########################
def create_data_model():
    """Stores the data for the problem"""
    data = {}
    # Locations in block unit
    _locations = \
            [(4, 4), # depot
             (2, 0), (8, 0), # locations to visit
             (0, 1), (1, 1),
             (5, 2), (7, 2),
             (3, 3), (6, 3),
             (5, 5), (8, 5),
             (1, 6), (2, 6),
             (3, 7), (6, 7),
             (0, 8), (7, 8)]
    _locations = mutate_locations(_locations)
    # Compute locations in meters using the block dimension defined as follow
    # Manhattan average block: 750ft x 264ft -> 228m x 80m
    # here we use: 114m x 80m city block
    # src: https://nyti.ms/2GDoRIe "NY Times: Know Your distance"
    data['locations'] = [(l[0] * 114, l[1] * 80) for l in _locations]
    data['num_locations'] = len(data['locations'])
    data['time_windows'] = \
          [(0, 0),
           (75, 85), (75, 85), # 1, 2
           (60, 70), (45, 55), # 3, 4
           (0, 8), (50, 60), # 5, 6
           (0, 10), (10, 20), # 7, 8
           (0, 10), (75, 85), # 9, 10
           (85, 95), (5, 15), # 11, 12
           (15, 25), (10, 20), # 13, 14
           (45, 55), (30, 40)] # 15, 16
    data['demands'] = \
          [0, # depot
           1, 1, # 1, 2
           2, 4, # 3, 4
           2, 4, # 5, 6
           8, 8, # 7, 8
           1, 2, # 9,10
           1, 2, # 11,12
           4, 4, # 13, 14
           8, 8] # 15, 16
    data['time_per_demand_unit'] = 5  # 5 minutes/unit
    data['num_vehicles'] = 6
    data['breaks'] = [(2, False), (2, False), (2, False), (2, False)]
    data['vehicle_capacity'] = 15
    data['vehicle_speed'] = 83  # Travel speed: 5km/h converted in m/min
    data['depot'] = 0
    return data


#######################
# Problem Constraints #
#######################
def manhattan_distance(position_1, position_2):
    """Computes the Manhattan distance between two points"""
    return (
        abs(position_1[0] - position_2[0]) + abs(position_1[1] - position_2[1]))


def create_distance_evaluator(data):
    """Creates callback to return distance between points."""
    _distances = {}
    # precompute distance between location to have distance callback in O(1)
    for from_node in range(data['num_locations']):
        _distances[from_node] = {}
        for to_node in range(data['num_locations']):
            if from_node == to_node:
                _distances[from_node][to_node] = 0
            else:
                _distances[from_node][to_node] = (manhattan_distance(
                    data['locations'][from_node], data['locations'][to_node]))

    def distance_evaluator(manager, from_node, to_node):
        """Returns the manhattan distance between the two nodes"""
        return _distances[manager.IndexToNode(from_node)][manager.IndexToNode(
            to_node)]

    return distance_evaluator


def create_demand_evaluator(data):
    """Creates callback to get demands at each location."""
    _demands = data['demands']

    def demand_evaluator(manager, node):
        """Returns the demand of the current node"""
        return _demands[manager.IndexToNode(node)]

    return demand_evaluator


def add_capacity_constraints(routing, data, demand_evaluator_index):
    """Adds capacity constraint"""
    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):
    """Creates callback to get total times between locations."""

    def service_time(data, node):
        """Gets the service time for the specified location."""
        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)))


###########
# Printer #
###########
def print_solution(data, manager, routing, assignment):  # pylint:disable=too-many-locals
    """Prints assignment on console"""
    print('Objective: {}'.format(assignment.ObjectiveValue()))

    print('Breaks:')
    intervals = assignment.IntervalVarContainer()
    for i in range(intervals.Size()):
        brk = intervals.Element(i)
        if brk.PerformedValue() == 1:
            print('{}: Start({}) Duration({})'.format(
                brk.Var().Name(),
                brk.StartValue(),
                brk.DurationValue()))
        else:
            print('{}: Unperformed'.format(brk.Var().Name()))

    total_distance = 0
    total_load = 0
    total_time = 0
    capacity_dimension = routing.GetDimensionOrDie('Capacity')
    time_dimension = routing.GetDimensionOrDie('Time')
    for vehicle_id in range(data['num_vehicles']):
        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))


########
# Main #
########
def main():
    """Entry point of the program"""
    # Instantiate the data problem.
    data = create_data_model()

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

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

    # Define weight of each edge
    distance_evaluator_index = routing.RegisterTransitCallback(
        partial(create_distance_evaluator(data), manager))
    routing.SetArcCostEvaluatorOfAllVehicles(distance_evaluator_index)

    # Add Capacity constraint
    demand_evaluator_index = routing.RegisterUnaryTransitCallback(
        partial(create_demand_evaluator(data), manager))
    add_capacity_constraints(routing, data, demand_evaluator_index)

    # Add Time Window constraint
    time_evaluator_index = routing.RegisterTransitCallback(
        partial(create_time_evaluator(data), manager))
    add_time_window_constraints(routing, manager, data, time_evaluator_index)

    # Add breaks
    time_dimension = routing.GetDimensionOrDie("Time")
    node_visit_transit = {}
    for n in range(routing.Size()):
        if n >= data['num_locations']:
            node_visit_transit[n] = 0
        else:
            node_visit_transit[n] = int(
                data['demands'][n] * data['time_per_demand_unit'])

    break_intervals = {}
    #for v in xrange(data['num_vehicles']):
    for v in [0]:
        vehicle_break = data['breaks'][v]
        break_intervals[v] = [
            routing.solver().FixedDurationIntervalVar(
                15, 100, vehicle_break[0], vehicle_break[1], 'Break for vehicle {}'.format(v))
            ]
        time_dimension.SetBreakIntervalsOfVehicle(
            break_intervals[v], v, node_visit_transit)

    # Setting first solution heuristic (cheapest addition).
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)  # pylint: disable=no-member
    # Solve the problem.
    assignment = routing.SolveWithParameters(search_parameters)
    print_solution(data, manager, routing, assignment)

    check_solution_output(data, manager, routing, assignment)



if __name__ == '__main__':
    main()

[Mizux]

1) Since you are working with time windows solver will return feasible range value which is harder to interpret since you can't just use all min/max of cumulVar and SlackVar.
i.e. the solver sometime can shift from one or two unit of time without changing the objective cost

see: https://github.com/google/or-tools/issues/1051#issuecomment-461033132

2) You can try to "force" the solver to return the earliest time instead of this "multi solution mixture" but I didn't see any AddVariableMinimizedByFinalizer in your code.

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

ref: https://developers.google.com/optimization/routing/vrptw#time_window_constraints

[ivarbrek]

1. In my code I use "Value" from Cumulvar and SlackVar and not Min/Max. Doesn't the solver ensure that the outputted "Values" satisfy slack(i) = cumul(j) - cumul(i) - transit(i, j)?
   After extending my assertion error output I can see that there exists combinations of time_i and slack_i values that satisfy the equation (e.g. time_i=8 and slack_i=21 or time_i=3 and slack_i=26).
   AssertionError: Time equation not in balance between nodes 5–4: 45!=3+16+21 indexes: 5-4
time_i range: (3, 8) transit_ij: 16 time_j range: (45, 45) slack_i range: (21, 26)

2. The mentioned problem still occurs when using AddVariableMinimizedByFinalizer

[KeremAslan]

I second that this line does not return the correct value:

slack_i = assignment.Value(time_dimension.SlackVar(i))

I have tested it with the following use case:

node | tw_start | tw_end
0    |      0   |  1000         
1    |    100   |    200

where the transit time from 0 -> 1 is 25.

Based on slack(i) = cumul(j) - cumul(i) - transit(i, j) I would expect a slack(0) to be 75. Yet assignment.value(slackVar) returns 0

Adding AddVariableMinimizedByFinalizer however seems to do the trick

[lperron]

Are you using the index manager in your check?

Le jeu. 8 oct. 2020 à 03:03, KeremAslan [e-mail] a écrit :

I second that this line does not return the correct value:

slack_i = assignment.Value(time_dimension.SlackVar(i))

I have tested it with the following use case:

node | tw_start | tw_end
0 | 0 | 1000
1 | 100 | 200

where the transit time from 0 -> 1 is 25.

Based on slack(i) = cumul(j) - cumul(i) - transit(i, j) I would expect a
slack(0) to be 75. Yet assignment.value(slackVar) returns 0

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub
https://github.com/google/or-tools/issues/2101#issuecomment-705271343,
or unsubscribe
https://github.com/notifications/unsubscribe-auth/ACUPL3KCMGA3KJKAJFMRMHTSJUFVXANCNFSM4OZV6ISQ
.

[KeremAslan]

Are you using the index manager in your check?

SlackVar I am accessing via index and CumulVar I am accessing via the index manager. That said in my test both the original index and indexManager.indexToNode(index) return the same value.

[lperron]

No, NodeToIndex(), you can only query the slack on indexes that come from
valid nodes.
Laurent Perron | Operations Research | [email protected] | (33) 1 42 68 53
00

Le jeu. 8 oct. 2020 à 08:36, KeremAslan [e-mail] a écrit :

Are you using the index manager in your check?

I am using it with the original index, but in my test both the original
index and indexManager.indexToNode(index) return the same value.

—
You are receiving this because you commented.
Reply to this email directly, view it on GitHub
https://github.com/google/or-tools/issues/2101#issuecomment-705363946,
or unsubscribe
https://github.com/notifications/unsubscribe-auth/ACUPL3PFOUXSRTZIUUVKSUDSJVMYRANCNFSM4OZV6ISQ
.

[KeremAslan]

I don't think I understand:

def print_solution(data, manager, routing, solution):
    """Prints solution on console."""
    time_dimension = routing.GetDimensionOrDie('Time')
    total_time = 0
    for vehicle_id in range(data['num_vehicles']):
        index = routing.Start(vehicle_id)
        plan_output = 'Route for vehicle {}:\n'.format(vehicle_id)
        while not routing.IsEnd(index):
            time_var = time_dimension.CumulVar(index)
            slack_var = time_dimension.SlackVar(manager.nodeToIndex(index))# <- this is what you mean?
            plan_output += '{0} Time({1},{2}) -> '.format(
                manager.IndexToNode(index), solution.Min(time_var),
                solution.Max(time_var))
            index = solution.Value(routing.NextVar(index))
        time_var = time_dimension.CumulVar(index)
        plan_output += '{0} Time({1},{2})\n'.format(manager.IndexToNode(index),
                                                    solution.Min(time_var),
                                                    solution.Max(time_var))
        plan_output += 'Time of the route: {}min\n'.format(
            solution.Min(time_var))
        print(plan_output)
        total_time += solution.Min(time_var)
    print('Total time of all routes: {}min'.format(total_time))

Are you saying that in this example I need to access slackVar via:

slack_var = time_dimension.SlackVar(manager.nodeToIndex(index))

and not

slack_var = time_dimension.SlackVar(manager.indexToNode(index))

?

[lperron]

No,

I meant calling dimension.GetSlackVar(1) may not be valid. You do not know
if index=1 point to a correct node.
Laurent Perron | Operations Research | [email protected] | (33) 1 42 68 53
00

Le jeu. 8 oct. 2020 à 09:06, KeremAslan [e-mail] a écrit :

I don't think I understand:

def print_solution(data, manager, routing, solution):
"""Prints solution on console."""
time_dimension = routing.GetDimensionOrDie('Time')
total_time = 0
for vehicle_id in range(data['num_vehicles']):
index = routing.Start(vehicle_id)
plan_output = 'Route for vehicle {}:n'.format(vehicle_id)
while not routing.IsEnd(index):
time_var = time_dimension.CumulVar(index)
slack_var = time_dimension.SlackVar(manager.nodeToIndex(index))# <- this is what you mean?
plan_output += '{0} Time({1},{2}) -> '.format(
manager.IndexToNode(index), solution.Min(time_var),
solution.Max(time_var))
index = solution.Value(routing.NextVar(index))
time_var = time_dimension.CumulVar(index)
plan_output += '{0} Time({1},{2})n'.format(manager.IndexToNode(index),
solution.Min(time_var),
solution.Max(time_var))
plan_output += 'Time of the route: {}minn'.format(
solution.Min(time_var))
print(plan_output)
total_time += solution.Min(time_var)
print('Total time of all routes: {}min'.format(total_time))

Are you saying that in this example I need to access slackVar via:

slack_var = time_dimension.SlackVar(manager.nodeToIndex(index))

and not

slack_var = time_dimension.SlackVar(manager.indexToNode(index))

?

—
You are receiving this because you commented.
Reply to this email directly, view it on GitHub
https://github.com/google/or-tools/issues/2101#issuecomment-705376726,
or unsubscribe
https://github.com/notifications/unsubscribe-auth/ACUPL3INVR3V2NSZDKVUDRTSJVQIZANCNFSM4OZV6ISQ
.