Routing - 3 features I do not know how to implement. Dimension based on another accumulated dimension, nodes that need another node to be visited earlier, and time intervals between node visits

[Marrit Schellekens]

I'm using google OR-tools in python to schedule a real-world routing problem. The problem is complex in that it has 2 different depots and 2 kinds of delivery goods. I already have a pretty reasonably working version with demands, maximum capacities, cost functions and time windows. There are however three different features I have no clue how to implement. So any tips or guidance to the right direction/examples/functions to look at would be appreciated. The three issues are:

1) Costs per 24 hours use

I would like to have a constant cost for every 24 hours that a vehicle is in use. I currently just have a constant cost for every vehicle that is in use. I basically want a dimension that is based on another accumulated dimension. I have seen functions for this in the C++ documentation, but couldn't find anything about in in python

2) Visit the correct depot before delivering goods.

Because I have two types of goods and two depots, I have to make sure to visit the correct depot before doing any deliveries of that delivery type. I figured I could do it with negative/positive demands, but that would mean a copy of the depot for each delivery, which would greatly increase the number of nodes. I was wondering if anyone had a more elegant suggestion.

3) Time intervals between node visits 

This is the most important one. I would like that some nodes have a certain time interval between visiting them (minimum and maximum). I could just choose a time for one of them, and base the time intervals based on that, but that would of course greatly limit the search space. So I was wondering if there was a way to specify the time interval between nodes instead of the time interval for one node.

[blind.lin...@gmail.com]

Marrit,

Taking point 3 first...

On Thu, Jul 16, 2020 at 12:14:13PM -0700, Marrit Schellekens wrote:
...
> *3) Time intervals between node visits *

> This is the most important one. I would like that some nodes have a certain
> time interval between visiting them (minimum and maximum). I could just
> choose a time for one of them, and base the time intervals based on that,
> but that would of course greatly limit the search space. So I was wondering
> if there was a way to specify the time interval between nodes instead of
> the time interval for one node.

I'm guessing a bit based on your description but,

Suppose there is a set of nodes A and another set of nodes B (the sets
might be identical) and you want to make sure that for every node in
set A, all nodes in set B are served either N time units earlier or
later. Further, suppose time is a dimension that has been set up
properly and is called "time".

My approach would be to just constrain the time. Make sure that the
time at A is either N units earlier or N units later than the time at
B. Then make sure that this constraint is applied if and only if the
same vehicle is serving A and B (assuming that is what you want).

This is just pseudo code, I haven't run it, but I think something like
this will work:

(in python)
```python

time_dimension = routing.GetDimensionOrDie("time")

for a_node in demand_subset_A:
a_index = manager.NodeToIndex(a_node)
for b_node in demand_subset_B:
if a_node == b_node:
# same node, so skip it
continue
b_index = manager.NodeToIndex(b_node)

# A comes before B by at least N
A_then_B_constraint = time_dimension.CumulVar(a_index) <=
time_dimension.CumulVar(b_index) - N

# or B comes before A by at least N
B_then_A_constraint = time_dimension.CumulVar(b_index) <=
time_dimension.CumulVar(a_index) - N

# one of the two constraints must hold
# if not using python, you have to combine using MakeSum and require equality using MakeEquality
# this also assumes that true is 1, false is 0,
# so true plus false and false plus true are both 1 (what we want),
# false plus false is zero (not allowed),
# true plus true is 2 (not allowed)

gap_constraint = A_then_B_constraint + B_then_A_constraint == 1

# conditional constraint...only apply if same vehicle
same_vehicle = routing.VehicleVar(a_index) == routing.VehicleVar(b_index)

conditional_constraint = solver.ConditionalExpression(
same_vehicle, gap_constraint, 1)

solver.Add(conditional_constraint>=1)

```

Caveat: I'm never confident about the ConditionalExpression semantics
until I run them and test out what is happening. The documentation is
super thin, and there are no official examples that I can find.

Hope that helps,
James

>
> --
> You received this message because you are subscribed to the Google Groups "or-tools-discuss" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to or-tools-discu...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/or-tools-discuss/27254864-d820-4d28-9ead-053e0b993b8do%40googlegroups.com.


--

James E. Marca
Activimetrics LLC

[blind.lin...@gmail.com]

Marrit,

Next taking point 2

On Thu, Jul 16, 2020 at 12:14:13PM -0700, Marrit Schellekens wrote:
...
>

> *2) Visit the correct depot before delivering goods. *

> Because I have two types of goods and two depots, I have to make sure to
> visit the correct depot before doing any deliveries of that delivery type.
> I figured I could do it with negative/positive demands, but that would mean
> a copy of the depot for each delivery, which would greatly increase the
> number of nodes. I was wondering if anyone had a more elegant suggestion.
>

So for this one I'm not sure what you are asking exactly. It seems
like you are wanting to set up a pickup and delivery problem. The PDP
problem *requires* one pickup per delivery. I've seen in the
documentation discussion of visiting a depot and reloading and then
delivering, but frankly I've never used that approach myself because
my problems always tend to be pickup and delivery.

So if you just use pickup and delivery, and do create a unique pickup
node (that is a dummy copy of the actual depot location) for every
delivery node, then problem solved.

But from your question, it appears that you do *not* want to go the
full PDP route, but would rather visit a depot, load up, and then
deliver to nodes until the vehicle is empty, but you want to make sure
that the correct depot is visited prior to visiting the nodes.

I'm really not sure about whether this will work in practice (I can
see issues with needing to reload, tracking dummy depot nodes, etc),
but my approach would be to use a "counting" dimension and make sure
that the value at the "correct" depot is less than the value at the
node. As with my suggestion for issue 3, you have to add a
conditional constraint making sure that the vehicles are the same,
because the count dimension is only meaningful for a single vehicle.

Assume for product A there is depot_a and nodes_a.

```python
count_dimension_name = 'count'
# count
routing.AddConstantDimension(
1, # increment by one every time
num_nodes, # max count is visit all the nodes
True, # set count to zero
count_dimension_name)
count_dimension = routing.GetDimensionOrDie(count_dimension_name)

# product A
# comes from some a depot. your implementation here is important
depot_a_idx = manager.NodeToIndex(depot_a_node)

# not sure about this; might be nonsensical if not careful about how
# depots are handled per vehicle
depot_count = count_dimension.CumulVar(depot_a_idx)

for a_node in nodes_a:
a_idx = manager.NodeToIndex(a_node)
a_count = count_dimension.CumulVar(a_idx)
# visit depot first
depot_first = solver.AddConstraint(depot_count < a_count)
# same vehicle
same_vehicle = routing.VehicleVar(depot_a_idx) == routing.VehicleVar(a_idx)
# require depot first only if same vehicle
conditional = solver.ConditionalExpression(same_vehicle, depot_first, 1)
solver.Add(conditional >= 1)

```

Again, if you just go with a full-blown pickup and delivery problem,
yes you will increase the numbers of nodes, but you also won't need
the above constraints, as the usual pickup and delivery constraints
will require:

same vehicle at pickup_depot_123 and delivery_node_123
time at pickup_depot_123 < time at delivery_node_123

Hope that helps,
James


--

[blind.lin...@gmail.com]

Marrit,

And finally, a suggestion for point 1

On Thu, Jul 16, 2020 at 12:14:13PM -0700, Marrit Schellekens wrote:
...
>

> *1) Costs per 24 hours use*

> I would like to have a constant cost for every 24 hours that a vehicle is
> in use. I currently just have a constant cost for every vehicle that is in
> use. I basically want a dimension that is based on another accumulated
> dimension. I have seen functions for this in the C++ documentation, but
> couldn't find anything about in in python
>

You can't use dependent dimensions in python, as you discovered.

I would suggest creating an explicit daily cost dimension. Regular
nodes add zero to the cost dimension, but create a dummy "end-of-day"
node, one per day per vehicle, that adds your daily cost. So handle
this logic in the callback function (check node type, if not an
end-of-day-node, return zero, otherwise return the daily cost).

Then add all the end of day nodes to the solver's minimization by
using a call to
routing.AddVariableMinimizedByFinalizer(daily_cost_dimension.CumulVar(day_i_vehicle_j_index))

Maybe that points you in the right direction?

Regards,

[Marrit Schellekens]

Thanks so much James for your time and detailed answers. I had a hard time finding documentation on the constraints method, so wasn't sure if it was what I need or if it was even supported in Python.  Sounds like that is indeed exactly what I need for 1. For 2, I started my modeling of doing pickup-and-delivery, but had a hard time getting it to work properly with my maximum capacities and demands, because I wanted to avoid doubling all my nodes with an additional depot 'dummy' node. However, now that I have a working version I'm very impressed with the low amount of time required for a decent route, so I'm not as scared of introducing extra nodes. I might switch over to doing full pickup-and-delivery. As for 3, that's a very smart little trick! I already have separate costs for different vehicle indexes, so this should be easy enough to implement!

I will be working on including these the coming few days, thank you so much for all your help! It is very much appreciated.

[blind.lin...@gmail.com]

Marrit,

You're welcome.

I was just thinking about point 1. I don't think you need to use an
extra cost dimension, but rather can just merge with whatever you're
already having the model minimize (distance or time).

What I mean is create the daily node as I described (one per vehicle
per day), allow them to be dropped (disjunction of zero) but require
them to be visited if the time dimension for a vehicle is greater than
a multiple of 24 hours. Then make the distance to the daily node
equal to the extra cost per day that you want to incur.

So for example, in the standard examples, you have this:

def create_data_model():
"""Stores the data for the problem."""
data = {}
data['distance_matrix'] = [
...
# set distance from all nodes to dummy end-of-day nodes as daily cost
]

...

# Instantiate the data problem.
data = create_data_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)
#
# this next line will correctly add "cost" of visiting end-of-day nodes
# if the data was set up correctly earlier
#
return data['distance_matrix'][from_node][to_node]


transit_callback_index = routing.RegisterTransitCallback(distance_callback)
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

> --
> You received this message because you are subscribed to the Google Groups "or-tools-discuss" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to or-tools-discu...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/or-tools-discuss/f9a6184a-61f3-4cf3-b8d3-23c5fd01369do%40googlegroups.com.

[Iván Babic]

Hi Marrit!

I'm struggling with a really similar case and I cant make point 2 work, I tried the solution as stated but the following error arose:

File "/.../python3.7/site-packages/ortools/constraint_solver/pywrapcp.py", line 453, in ConditionalExpression
    return _pywrapcp.Solver_ConditionalExpression(self, condition, expr, unperformed_value)
SystemError: <built-in function Solver_ConditionalExpression> returned NULL without setting an error

I figured out that the condition depot_first should be added to the solver but declared as:

depot_first = count_dimension.CumulVar(depot_type_idx) <= count_dimension.CumulVar(start_index)




But now the solver doesn't find a solution when I add the condition

same_vehicle = routing.VehicleVar(depot_type_idx) == routing.VehicleVar(start_index)
depot_first = time_dimension.CumulVar(depot_type_idx) <= time_dimension.CumulVar(start_index)
conditional = routing.solver().ConditionalExpression(same_vehicle, depot_first, 1)
routing.solver().Add(conditional >= 1)




Is there a way of debugging what is happening? I will post my code here, which doesn't find solutions unless the conditional is  set to be  >=0




"""Capacited Vehicles Routing Problem (CVRP)."""

from __future__ import print_function
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp

from collections import namedtuple

def create_data_model():
    """Stores the data for the problem."""
    data = {}

    data['time_matrix'] = [
        [0, 1, 1, 2, 3, 6, 2, 3, 4],
        [1, 0, 1, 3, 2, 6, 8, 4, 3],
        [1, 1, 0, 11, 10, 6, 3, 9, 1],
        [8, 3, 11, 0, 1, 7, 10, 6, 2],
        [7, 2, 10, 1, 0, 6, 9, 4, 3],
        [3, 6, 6, 7, 6, 0, 2, 3, 3],
        [6, 8, 3, 10, 9, 2, 0, 6, 6],
        [2, 4, 9, 6, 4, 3, 6, 0, 2],
        [2, 4, 9, 6, 4, 3, 6, 2, 0]
    ]
    PickUpDeliveryInfo = namedtuple('PickUpDeliveryInfo', ['start_node', 'end_node', 'load_type', 'duration'])
    pickups_deliveries_test = [
        [3, 4, 0, 1],
        [5, 6, 0, 1]
    ]
    data['pickups_deliveries'] = [
        PickUpDeliveryInfo(start_node, end_node, load_type, duration)
        for start_node, end_node, load_type, duration in pickups_deliveries_test
    ]

    data['type_depots'] = {
        0: 1,
        1: 2
    }
    data['max_time'] = 100
    data['capacity'] = 100
    data['demands'] = {
        'TypeA': [x * data['capacity']/10 for x in [0,  0,  0,  10, -10, 10, -10,  0,  0]],
        'TypeB': [x * data['capacity']/10 for x in [0,  0,  0,  0,  0,  0,  0,  0, 0]]
    }

    data['num_vehicles'] = 1
    data['depot'] = 0
    return data


def get_type_depot(data, load_type):
    return data['type_depots'][load_type]


def print_node_dimension(routing, solution, index, dimension, end=False):
    cumul_var = dimension.CumulVar(index)
    plan_output = f'Cumul[{solution.Min(cumul_var)},{solution.Max(cumul_var)}] {solution.Value(cumul_var)}'
    if not routing.IsStart(index) and not end:
        trans_var = dimension.TransitVar(index)
        plan_output += f'; Transit[{solution.Min(trans_var)},{solution.Max(trans_var)}] {solution.Value(trans_var)}'
        slack_var = dimension.SlackVar(index)
        plan_output += f'; Slack[{solution.Min(slack_var)},{solution.Max(slack_var)}] {solution.Value(slack_var)}'
    print(plan_output)


def print_solution(data, manager, routing, solution):
    """Prints solution on console."""
    total_distance = 0
    total_load = 0
    capacityA: pywrapcp.RoutingDimension = routing.GetDimensionOrDie('CapacityA')
    capacityB: pywrapcp.RoutingDimension = routing.GetDimensionOrDie('CapacityB')
    count: pywrapcp.RoutingDimension = routing.GetDimensionOrDie('count')
    for vehicle_id in range(data['num_vehicles']):
        index = routing.Start(vehicle_id)
        print('Route for vehicle {}:\n'.format(vehicle_id))
        route_distance = 0
        route_load_A = 0
        while not routing.IsEnd(index):
            print(f'[[[[{manager.IndexToNode(index)}]]]]: ')
            print('CapacityA')
            print_node_dimension(routing, solution, index, capacityA)
            print('CapacityB')
            print_node_dimension(routing, solution, index, capacityB)
            print('Count')
            print_node_dimension(routing, solution, index, count, True)
            print(' -> ')
            index = solution.Value(routing.NextVar(index))
        print(f'[[[[{manager.IndexToNode(index)}]]]]: ')
        print_node_dimension(routing, solution, index, capacityA, True)
        print_node_dimension(routing, solution, index, capacityB, True)
        total_distance += route_distance
        total_load += route_load_A
    print('Total distance of all routes: {}m'.format(total_distance))
    print('Total load of all routes: {}'.format(total_load))


def main():
    """Solve the CVRP problem."""

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

    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)
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

    # Add Time constraint.
    time = 'Time'
    routing.AddDimension(
        transit_callback_index,
        data['max_time'],  # allow waiting time
        data['max_time'],  # maximum time per vehicle
        True,  # Don't force start cumul to zero.
        time)
    time_dimension = routing.GetDimensionOrDie(time)

    count_dimension_name = 'count'
    # count
    routing.AddConstantDimension(

        1,# increment by one every time
        len(data['time_matrix']),  # max count is visit all the nodes

        True,  # set count to zero
        count_dimension_name)
    count_dimension = routing.GetDimensionOrDie(count_dimension_name)

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

    # Add Capacity constraint.
    def a_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']['TypeA'][from_node]

    a_demand_callback_index = routing.RegisterUnaryTransitCallback(a_demand_callback)
    routing.AddDimension(
        a_demand_callback_index,
        0,
        data['capacity'],  # vehicle maximum capacities
        True,  # start cumul to zero
        'CapacityA')

    # Add Capacity constraint.
    def b_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']['TypeB'][from_node]

    b_demand_callback_index = routing.RegisterUnaryTransitCallback(b_demand_callback)
    routing.AddDimension(
        b_demand_callback_index,
        0,  # null capacity slack
        data['capacity'],  # vehicle maximum capacities
        True,  # start cumul to zero
        'CapacityB')

    a_dimension = routing.GetDimensionOrDie('CapacityA')
    b_dimension = routing.GetDimensionOrDie('CapacityB')

    # for node_index in range(len(data['time_matrix'][0])):
    #     index = manager.NodeToIndex(node_index)
    #     if node_index == 1:
    #         a_dimension.SlackVar(index).SetRange(0, data['capacity'])
    #         b_dimension.SlackVar(index).SetValue(0)
    #     elif node_index == 2:
    #         a_dimension.SlackVar(index).SetValue(0)
    #         b_dimension.SlackVar(index).SetRange(0, data['capacity'])
    #     elif node_index != 0:
    #         a_dimension.SlackVar(index).SetValue(0)
    #         b_dimension.SlackVar(index).SetValue(0)

    for pd_info in data['pickups_deliveries']:
        start_index = manager.NodeToIndex(pd_info.start_node)
        end_index = manager.NodeToIndex(pd_info.end_node)

        routing.solver().Add(routing.NextVar(start_index) == end_index)
        routing.AddPickupAndDelivery(start_index, end_index)
        routing.solver().Add(
            routing.VehicleVar(start_index) == routing.VehicleVar(
                end_index))
        routing.solver().Add(
            time_dimension.CumulVar(start_index) <=
            time_dimension.CumulVar(end_index))
        depot_type_idx = manager.NodeToIndex(get_type_depot(data, pd_info.load_type))
        same_vehicle = routing.VehicleVar(depot_type_idx) == routing.VehicleVar(start_index)
        depot_first = time_dimension.CumulVar(depot_type_idx) <= time_dimension.CumulVar(start_index)
        conditional = routing.solver().ConditionalExpression(same_vehicle, depot_first, 1)
        routing.solver().Add(conditional >= 1)

    for location_index in range(routing.nodes()):
        index = manager.NodeToIndex(location_index)
        if routing.IsEnd(index) or routing.IsStart(index):
            continue
        routing.AddToAssignment(a_dimension.TransitVar(index))
        routing.AddToAssignment(a_dimension.SlackVar(index))
        routing.AddToAssignment(b_dimension.TransitVar(index))
        routing.AddToAssignment(b_dimension.SlackVar(index))

    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC

    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PARALLEL_CHEAPEST_INSERTION)

    search_parameters.time_limit.seconds = 5  # Amount of time you want to give the solver, for searching for better solutions
    search_parameters.local_search_metaheuristic = (
        routing_enums_pb2.LocalSearchMetaheuristic.SIMULATED_ANNEALING
    )
    search_parameters.log_search = True
    # Solve the problem.
    solution = routing.SolveWithParameters(search_parameters)
    print(solution)
    # Print solution on console.
    if solution:
        print_solution(data, manager, routing, solution)


if __name__ == '__main__':
    main()

[Marrit Schellekens]

I have not yet tried implementing point 2, when I do I'll reply here.

[Marrit Schellekens]

Replying here in case anyone ever has the same question. I managed to implement this successfully. My final code looked like this:

    for set in data['sets']:

        #For all possible combinations in a set:

        for i,index in enumerate(set):

            for j in range((i+1), len(set)):

                # A comes before B by at least N

                A_then_B_constraint = time_dimension.CumulVar(set[i]) <= time_dimension.CumulVar(set[j]) - constants['MIN_BETWEEN_ORDERS']

                # or B comes before A by at least N

                B_then_A_constraint = time_dimension.CumulVar(set[j]) <= time_dimension.CumulVar(set[i]) - constants['MIN_BETWEEN_ORDERS']

                gap_constraint = A_then_B_constraint + B_then_A_constraint == 1

                routing.solver().Add(gap_constraint)

> To unsubscribe from this group and stop receiving emails from it, send an email to or-tools...@googlegroups.com.