OR-tools routing fails to find 'trivial' optimal solution in the presence of order-constraints
<p>In order to better understand the constraint-programming behind or-tools routing, I created a toy example of a depot and 4 other nodes configured to allow for two routes.</p>

<p><a href="http://i.stack.imgur.com/SiZj1.jpg" rel="nofollow"><img src="http://i.stack.imgur.com/SiZj1.jpg" alt="enter image description here"></a></p>

<p>The idea is that a vehicle travels from the depot <code>0</code> to <code>1</code>, then either picks <code>2</code> or <code>3</code>, continues to <code>4</code> and returns to the depot <code>0</code>; the vehicle either picks the green or red path. My actual problem is more complicated and has multiple vehicles, but has similar constraints.</p>

<p>For this example I created a euclidean-distance function for the cost: </p>



<pre class="lang-py prettyprint-override"><code>class Distances:
    def __init__(self):
        self.locations = [
            [-1,  0], # source
            [ 0, -1], # waypoint 1
            [ 0,  1], # waypoint 2
            [ 1,  0], # destination
        ]

    def __len__(self):
        return len(self.locations) + 1

    def dist(self, x, y):
        return int(10000 * math.sqrt(
            (x[0] - y[0]) ** 2 +
            (x[1] - y[1]) ** 2))

    def __call__(self, i, j):
        if i == 0 and j == 0:
            return 0
        if j == 0 or i == 0:
            return 1 # very small distance between depot and non-depot, simulating 0
        return self.dist(self.locations[i - 1], self.locations[j - 1])


distance = Distances()
</code></pre>

<p>And a l0 distance function to constraint the order:</p>

<pre class="lang-py prettyprint-override"><code># l0-distance to add order constraints
class Order:
    def __call__(self, i, j):
        return 0 if i == j else 1


order = Order()
</code></pre>

<p>Then I create the model and try to solve this:</p>

<pre class="lang-py prettyprint-override"><code>search_parameters = pywrapcp.RoutingModel.DefaultSearchParameters()
search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.ALL_UNPERFORMED)

search_parameters.local_search_metaheuristic = routing_enums_pb2.LocalSearchMetaheuristic.SIMULATED_ANNEALING
search_parameters.time_limit_ms = 3000

routing = pywrapcp.RoutingModel(len(distance), 1)

routing.SetArcCostEvaluatorOfAllVehicles(distance)
routing.SetDepot(0)
solver = routing.solver()

routing.AddDimension(order, int(1e18), int(1e18), True, "order")

# Since `ALL_UNPERFORMED` is used, each node must be allowed inactive
order_dimension = routing.GetDimensionOrDie("order")
routing.AddDisjunction([1], int(1e10))
routing.AddDisjunction([2, 3], int(1e10))
routing.AddDisjunction([4], int(1e10))

solver.AddConstraint(order_dimension.CumulVar(1) &lt;= order_dimension.CumulVar(2))
solver.AddConstraint(order_dimension.CumulVar(1) &lt;= order_dimension.CumulVar(3))

solver.AddConstraint(order_dimension.CumulVar(2) &lt;= order_dimension.CumulVar(4))
solver.AddConstraint(order_dimension.CumulVar(3) &lt;= order_dimension.CumulVar(4))

# routing.AddPickupAndDelivery(1, 2)
# routing.AddPickupAndDelivery(1, 3)
# routing.AddPickupAndDelivery(2, 4)
# routing.AddPickupAndDelivery(3, 4)

routing.CloseModelWithParameters(search_parameters)
assignment = routing.SolveWithParameters(search_parameters)

if assignment is not None:
    print('cost: ' + str(assignment.ObjectiveValue()))
    route = []
    index = routing.Start(0)
    while not routing.IsEnd(index):
        route.append(routing.IndexToNode(index))
        index = assignment.Value(routing.NextVar(index))
    for node in route:
        print(' - {:2d}'.format(node))
else:
    print('nothing found')
</code></pre>

<p>So <code>[1]</code> and <code>[4]</code> are disjunctions to allow for <code>ALL_UNPERFORMED</code> first solution to work, and the disjunction <code>[2, 3]</code> to state that either the green or red path should be chosen.</p>

<p>With these disjunctions the solver finds a solution, but if I add that <code>2</code>  and <code>3</code> should be visited after <code>1</code> and before <code>4</code>, the solver does not visit <code>2</code> or <code>3</code> at all. Why is this the case? Why can't the solver find the more optimal route <code>0 -&gt; 1 -&gt; 2/3 -&gt; 4 -&gt; 0</code> avoiding the <code>int(1e10)</code> disjunction penalty for <code>[2,3]</code>?</p>

<p><strong>EDIT:</strong></p>

<p>Soft-constraining the order-constraints by removing them and adding to <code>Distance.__call__</code>:</p>

<pre class="lang-py prettyprint-override"><code>if (i == 2 or j == 1) or (i == 3 or j == 1) or (i == 4 or j == 2) or (i == 4 or j == 3):
    return int(1e10)
</code></pre>

<p>to penalize a disallowed order, results in the route <code>0 -&gt; 2 -&gt; 1 -&gt; 4 -&gt; 0</code>. So I wonder why or-tools won't swap <code>1</code> and <code>2</code>, even when explicitly enabling <code>use_swap_active</code> and <code>use_relocate_neighbors</code> in <code>search_parameters.local_search_operators</code>.</p>

<p><strong>NOTE:</strong> Failed because it should have been:</p>

<pre class="lang-py prettyprint-override"><code>if (i == 2 and j == 1) or (i == 3 and j == 1) or (i == 4 and j == 2) or (i == 4 and j == 3):
    return int(1e10)
</code></pre>

<p>Concluding: the search-space is small, a better solution is in the neighborhood of <code>use_relocate_neighbors</code> of the returned solution, yet or-tools does not find it. Why?</p>

<h1>All code:</h1>

<pre class="lang-py prettyprint-override"><code>import pandas
import os.path

import numpy
import math
from ortools.constraint_solver import pywrapcp
from ortools.constraint_solver import routing_enums_pb2


class Distances:
    def __init__(self):
        self.locations = [
            [-1,  0], # source
            [ 0, -1], # waypoint 1
            [ 0,  1], # waypoint 2
            [ 1,  0], # destination
        ]

    def __len__(self):
        return len(self.locations) + 1

    def dist(self, x, y):
        return int(10000 * math.sqrt(
            (x[0] - y[0]) ** 2 +
            (x[1] - y[1]) ** 2))

    def __call__(self, i, j):
        if i == 0 and j == 0:
            return 0
        if j == 0 or i == 0:
            return 1 # very small distance between depot and non-depot, simulating 0
        return self.dist(self.locations[i - 1], self.locations[j - 1])


distance = Distances()


# l0-distance to add order constraints
class Order:
    def __call__(self, i, j):
        return 0 if i == j else 1


order = Order()

search_parameters = pywrapcp.RoutingModel.DefaultSearchParameters()
search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.ALL_UNPERFORMED)

search_parameters.local_search_metaheuristic = routing_enums_pb2.LocalSearchMetaheuristic.SIMULATED_ANNEALING
search_parameters.time_limit_ms = 3000

routing = pywrapcp.RoutingModel(len(distance), 1)

routing.SetArcCostEvaluatorOfAllVehicles(distance)
routing.SetDepot(0)
solver = routing.solver()

routing.AddDimension(order, int(1e18), int(1e18), True, "order")

# Since `ALL_UNPERFORMED` is used, each node must be allowed inactive
order_dimension = routing.GetDimensionOrDie("order")
routing.AddDisjunction([1], int(1e10))
routing.AddDisjunction([2, 3], int(1e10))
routing.AddDisjunction([4], int(1e10))

solver.AddConstraint(
    (routing.ActiveVar(2) == 0)
    or
    (order_dimension.CumulVar(1) &lt;= order_dimension.CumulVar(2))
)
solver.AddConstraint(
    (routing.ActiveVar(3) == 0)
    or
    (order_dimension.CumulVar(1) &lt;= order_dimension.CumulVar(3))
)


solver.AddConstraint(
    (routing.ActiveVar(2) == 0)
    or
    (order_dimension.CumulVar(2) &lt;= order_dimension.CumulVar(4))
)
solver.AddConstraint(
    (routing.ActiveVar(3) == 0)
    or
    (order_dimension.CumulVar(3) &lt;= order_dimension.CumulVar(4))
)

# routing.AddPickupAndDelivery(1, 2)
# routing.AddPickupAndDelivery(1, 3)
# routing.AddPickupAndDelivery(2, 4)
# routing.AddPickupAndDelivery(3, 4)

routing.CloseModelWithParameters(search_parameters)
assignment = routing.SolveWithParameters(search_parameters)

if assignment is not None:
    print('cost: ' + str(assignment.ObjectiveValue()))
    route = []
    index = routing.Start(0)
    while not routing.IsEnd(index):
        route.append(routing.IndexToNode(index))
        index = assignment.Value(routing.NextVar(index))
    for node in route:
        print(' - {:2d}'.format(node))
else:
    print('nothing found')
</code></pre>

<p><strong>@furnon</strong> On github answered my question via the github-issues list: <a href="https://github.com/google/or-tools/issues/252#issuecomment-249646587" rel="nofollow">https://github.com/google/or-tools/issues/252#issuecomment-249646587</a></p>



<p>First, constraint programming performs better on tighter constraints, I guess some things are searched for exaustively. In particular, I had to limit the order-dimension:</p>

<pre class="lang-py prettyprint-override"><code>routing.AddDimension(order, int(1e18), int(1e18), True, "order")
</code></pre>

<p>is better constrained via</p>

<pre class="lang-py prettyprint-override"><code>routing.AddDimension(order, len(distance) + 1 ,len(distance) + 1, True, "order")
</code></pre>

<p>Subsequently, the check on whether <code>2</code> or <code>3</code> is active is not needed, hence we can simplify the order-constrains like this:</p>

<pre class="lang-py prettyprint-override"><code>solver.AddConstraint(order_dimension.CumulVar(1) &lt;= order_dimension.CumulVar(2))
solver.AddConstraint(order_dimension.CumulVar(1) &lt;= order_dimension.CumulVar(3))
solver.AddConstraint(order_dimension.CumulVar(2) &lt;= order_dimension.CumulVar(4))
solver.AddConstraint(order_dimension.CumulVar(3) &lt;= order_dimension.CumulVar(4))
</code></pre>

<p>As I already did in the inline version, but not in the all-code version. Now a feasible solution is returned.</p>

python or-tools