# Subtour elimination constraints

**URL:** <https://discourse.julialang.org/t/subtour-elimination-constraints/106902>\
**Category:** Optimization (Mathematical)\
**Tags:** jump, cplex, gurobi\
**Created:** [November 29, 2023, 2:29pm UTC](https://discourse.julialang.org/t/subtour-elimination-constraints/106902 "2023-11-29T14:29:29Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![Victor\_Neto](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/victor_neto/32/204390_2.png) [@Victor\_Neto](https://discourse.julialang.org/u/Victor_Neto)\
**Post date:** [November 29, 2023, 2:29pm UTC](https://discourse.julialang.org/t/subtour-elimination-constraints/106902/1 "2023-11-29T14:29:30Z")

</div>

Dear colleagues,  
I am studying the school bus routing problem: [Redirecting](https://doi.org/10.1016/j.ejor.2013.02.025).  
I am reproducing this mixed-integer linear programming formulation. However, I was replacing the sub-tour constraints (3) with a flow constraint (as proposed by Miller, Tucker, and Zemlin, 1960).  
I’ve modeled this constraint as follows:

@constraint(model\_sbrp, [i in V, j in V, k in 1:n], u[i] - u[j] + n\*x[i,j,k] \>= n - 1)

Nevertheless, an error message was returned (out of bounds).  
I need some assistance to tackle this problem. Could someone help me, please?

```julia
@variable(model_sbrp, x[i in V, j in V, k = 1:n, i != j], Bin) 
@variable(model_sbrp, y[i in V, k = 1:n], Bin)
@variable(model_sbrp, z[i in V, l in S, k = 1:n], Bin) 
@variable(model_sbrp, u[i in V] >= 0 ) 

@objective(model_sbrp, Min, sum(c[i,j] * x[i,j,k] for i in V, j in V, k in 1:n if i != j))

for i in V
    for k in 1:n
        @constraint(model_sbrp, sum(x[i, j, k] for j in V if i!=j) == sum(x[j, i, k] for j in V if i!=j))
        @constraint(model_sbrp, sum(x[i, j, k] for j in V if i!=j) == y[i, k])
    end
end

@constraint(model_sbrp, [i in V, i != 1],sum(y[i,k] for k = 1:n) <= 1)
@constraint(model_sbrp, [l in S, i in V], sum(z[i,l,k] for k = 1:n) <= s[i,l])
@constraint(model_sbrp, [k in 1:n], sum(z[i,l,k] for i in V, l in S) <= C)
@constraint(model_sbrp, [i in V, l in S, k in 1:n], z[i,l,k] <= y[i,k])
@constraint(model_sbrp, [l in S], sum(z[i,l,k] for i in V, k = 1:n) == 1)
@constraint(model_sbrp, [i in V, j in V, k in 1:n], u[i] - u[j] + n*x[i,j,k] >= n - 1)
@constraint(model_sbrp, u[1] == 1)

```

i’m sending a small-sized test instance

C = 10

q = 8

c=  
[0.0 3.16 1.41 1.0 2.24;  
3.16 0.0 4.0 3.61 4.12;  
1.41 4.0 0.0 2.24 1.0;  
1.0 3.61 2.24 0.0 3.16;  
2.24 4.12 1.0 3.16 0.0]  
s=  
[0 0 0 0 0 0 0 0 0 0;  
0 0 1 1 1 0 0 0 1 0;  
0 0 0 1 0 0 1 1 0 1;  
1 0 1 0 0 0 0 1 1 0;  
0 1 0 0 0 1 1 1 0 0]  
i = 1  
V = [1 2 3 4 5]  
S = size(s,2)  
n = length(V)

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [November 29, 2023, 3:01pm UTC](https://discourse.julialang.org/t/subtour-elimination-constraints/106902/2 "2023-11-29T15:01:19Z")

</div>

> [@Victor\_Neto](#):
>
> `u[i] - u[j] + n*x[i,j,k] >= n - 1`

Suppose `x[i,j,k]` is zero, then `u[i] - u[j] >= n-1` i.e. `u[i] >= u[j]+n-1`. This inequality seems to go in the wrong direcction, as we would want the `u[...]`s to be close together.

---

<div class="post-metadata">

**Author:** ![Victor\_Neto](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/victor_neto/32/204390_2.png) [@Victor\_Neto](https://discourse.julialang.org/u/Victor_Neto)\
**Post date:** [November 29, 2023, 5:20pm UTC](https://discourse.julialang.org/t/subtour-elimination-constraints/106902/3 "2023-11-29T17:20:56Z")

</div>

Dear Dan,  
yes, you’re correct, but when i change for this constraints

u[i] - u[j] + n\*x[i,j,k] \<= n - 1

the model still does not eliminate the subtours.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [November 29, 2023, 5:51pm UTC](https://discourse.julialang.org/t/subtour-elimination-constraints/106902/4 "2023-11-29T17:51:35Z")

</div>

`u[...]` depends just on the vertex in the bus tour and not on the bus number. In order to ensure buses are not constrained to go all on the same path, each bus needs its own `u[...]` sequence (and associated constraints).

---

<div class="post-metadata">

**Author:** ![Victor\_Neto](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/victor_neto/32/204390_2.png) [@Victor\_Neto](https://discourse.julialang.org/u/Victor_Neto)\
**Post date:** [November 29, 2023, 6:47pm UTC](https://discourse.julialang.org/t/subtour-elimination-constraints/106902/5 "2023-11-29T18:47:16Z")

</div>

Dear Dan  
I’m working on a MIP formulation presented in Section 3 on the paper sent as an attachment and I’m just changing constraint n 3 for Miller–Tucker–Zemlin constraint. This is proposed by the paper itself on topic 5.2 page 7, no need to include constraint.  
Thanks for your help and attention.

(Attachment Original\_A metaheuristic for the school bus routing problem with bus stop selection.pdf is missing)

---

<div class="post-metadata">

**Author:** ![Victor\_Neto](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/victor_neto/32/204390_2.png) [@Victor\_Neto](https://discourse.julialang.org/u/Victor_Neto)\
**Post date:** [November 29, 2023, 6:51pm UTC](https://discourse.julialang.org/t/subtour-elimination-constraints/106902/6 "2023-11-29T18:51:57Z")

</div>

[https://www.sciencedirect.com/science/article/abs/pii/S0377221713001586](https://www.sciencedirect.com/science/article/abs/pii/S0377221713001586)

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [November 29, 2023, 6:55pm UTC](https://discourse.julialang.org/t/subtour-elimination-constraints/106902/7 "2023-11-29T18:55:33Z")

</div>

The constraint in the papar, has a ‘for all k’, which means it is actually k independent constraints, one for each bus (to eliminate subtours for the bus). The use of `u[...]` to achieve this should again have `k` independent `u[...]`s, one for each bus.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [November 29, 2023, 7:20pm UTC](https://discourse.julialang.org/t/subtour-elimination-constraints/106902/8 "2023-11-29T19:20:10Z")

</div>

Hi @Victor_Neto, just a comment: you can make things easier by providing a reproducible example. This should be a single block of code that someone can copy-paste and reproduce your issue.

Here’s how I would write your current model.

```julia
using JuMP, Gurobi
c = [
    0.0 3.16 1.41 1.0 2.24
    3.16 0.0 4.0 3.61 4.12
    1.41 4.0 0.0 2.24 1.0
    1.0 3.61 2.24 0.0 3.16
    2.24 4.12 1.0 3.16 0.0
]
s = [
    0 0 0 0 0 0 0 0 0 0
    0 0 1 1 1 0 0 0 1 0
    0 0 0 1 0 0 1 1 0 1
    1 0 1 0 0 0 0 1 1 0
    0 1 0 0 0 1 1 1 0 0
]
V = [1, 2, 3, 4, 5]
C, S, n = 10, size(s, 2), length(V)
model = Model(Gurobi.Optimizer)
@variables(model, begin
    x[i in V, j in V, k in 1:n; i != j], Bin
    y[i in V, k in 1:n], Bin
    z[i in V, l in S, k in 1:n], Bin
    u[i in V] >= 0
end) 
@objective(
    model,
    Min,
    sum(c[i, j] * x[i, j, k] for i in V, j in V, k in 1:n if i != j),
)
@constraints(model, begin
    [i in V, k in 1:n], sum(x[i, j, k] for j in V if i!=j) == sum(x[j, i, k] for j in V if i!=j)
    [i in V, k in 1:n], sum(x[i, j, k] for j in V if i!=j) == y[i, k]
    [i in V; i != 1], sum(y[i, :]) <= 1
    [i in V, l in S], sum(z[i, l, :]) <= s[i, l]
    [k in 1:n], sum(z[:, :, k]) <= C
    [i in V, l in S, k in 1:n], z[i, l, k] <= y[i, k]
    [l in S], sum(z[:, l, :]) == 1
    [i in V, j in V, k in 1:n], u[i] - u[j] + n * x[i, j, k] >= n - 1
    u[1] == 1
end)
optimize!(model)

```

But agree with @Dan, you need a few tweaks.

---

<div class="post-metadata">

**Author:** ![Victor\_Neto](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/victor_neto/32/204390_2.png) [@Victor\_Neto](https://discourse.julialang.org/u/Victor_Neto)\
**Post date:** [November 29, 2023, 7:43pm UTC](https://discourse.julialang.org/t/subtour-elimination-constraints/106902/9 "2023-11-29T19:43:44Z")

</div>

Thanks for the tip Oscar, I’m new here.  
I’ve been working on it since last week and I don’t know how to fix this code. I tried different ways but always found an error.

---

<div class="post-metadata">

**Author:** ![Victor\_Neto](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/victor_neto/32/204390_2.png) [@Victor\_Neto](https://discourse.julialang.org/u/Victor_Neto)\
**Post date:** [November 29, 2023, 7:52pm UTC](https://discourse.julialang.org/t/subtour-elimination-constraints/106902/10 "2023-11-29T19:52:10Z")

</div>

Is this correct ? @variable(model\_sbrp, u[i in V] \>= 0 )

V is a set of potentials stops

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [November 29, 2023, 8:09pm UTC](https://discourse.julialang.org/t/subtour-elimination-constraints/106902/11 "2023-11-29T20:09:35Z")

</div>

I don’t understand what the `y` and `z` variables are doing, but this should point you in the right direction. (Assuming `k` is the different busses, that each must do a tour?)

```Julia
using JuMP, Gurobi
c = [
    0.0 3.16 1.41 1.0 2.24
    3.16 0.0 4.0 3.61 4.12
    1.41 4.0 0.0 2.24 1.0
    1.0 3.61 2.24 0.0 3.16
    2.24 4.12 1.0 3.16 0.0
]
s = [
    0 0 0 0 0 0 0 0 0 0
    0 0 1 1 1 0 0 0 1 0
    0 0 0 1 0 0 1 1 0 1
    1 0 1 0 0 0 0 1 1 0
    0 1 0 0 0 1 1 1 0 0
]
V = 1:size(c, 1)
C, S, n = 10, 1:size(s, 2), length(V)
model = Model(Gurobi.Optimizer)
@variables(model, begin
    x[i in V, j in V, k in 1:n], Bin
    y[i in V, k in 1:n], Bin
    z[i in V, l in S, k in 1:n], Bin
    2 <= u[i in V, k in 1:n] <= length(V)
end) 
fix.(u[1, :], 1; force = true)
@objective(
    model,
    Min,
    sum(c[i, j] * x[i, j, k] for i in V, j in V, k in 1:n),
)
@constraints(model, begin
    # For all nodes and routes, flow out == flow in
    [i in V, k in 1:n], sum(x[i, :, k]) == sum(x[:, i, k])
    # Flow in and out of other nodes is 1
    [i in V; i != 1], sum(x[i, :, :]) == 1
    [i in V; i != 1], sum(x[:, i, :]) == 1
    # MTZ constraints
    [i in V, j in V, k in 1:n; i != 1 && j != 1],
        u[i, k] - u[j, k] + 1 <= (n - 1) * (1 - x[i, j, k])
    # y[i, k] is 1 if node i is in route k
    [i in V, k in 1:n], sum(x[i, :, k]) == y[i, k]
    # each node can be in at most one route, except node 1, which is in every 
    # route
    [i in V; i != 1], sum(y[i, :]) <= 1
    # I don't understand these constraints
    [i in V, l in S], sum(z[i, l, :]) <= s[i, l]
    [k in 1:n], sum(z[:, :, k]) <= C
    [i in V, l in S, k in 1:n], z[i, l, k] <= y[i, k]
    [l in S], sum(z[:, l, :]) == 1
end)
optimize!(model)
for k in 1:n
    print("Route $k: 1")
    i = 1
    while true
        for j in V
            if value(x[i, j, k]) > 0.5
                print(" => $j")
                i = j
                break
            end
        end
        if i == 1
            println()
            break
        end
    end
end

```

---

<div class="post-metadata">

**Author:** ![Victor\_Neto](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/victor_neto/32/204390_2.png) [@Victor\_Neto](https://discourse.julialang.org/u/Victor_Neto)\
**Post date:** [November 29, 2023, 9:08pm UTC](https://discourse.julialang.org/t/subtour-elimination-constraints/106902/12 "2023-11-29T21:08:29Z")

</div>

variable y : 1 if bus k visits stop i, 0 otherwise  
variable z: 1 if student l is picked up by bus k at stop i, 0 otherwise

about the bus

for each k a route must start from index 1 (school) and also end at the school.

The problem I’m having is that one route is being done correctly (starting and finishing at school) but the others are not. generating the subroutes

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [November 29, 2023, 9:33pm UTC](https://discourse.julialang.org/t/subtour-elimination-constraints/106902/13 "2023-11-29T21:33:46Z")

</div>

Running odow’s code, I get:

```julia
Route 1: 1 => 1
Route 2: 1 => 1
Route 3: 1 => 1
Route 4: 1 => 1
Route 5: 1 => 3 => 5 => 2 => 4 => 1

```

What did you get?

(in some sense this is optimal, as Bus 5 takes all the students and the rest of the buses wait. The objective is to minimize travel costs while taking all students).

ADDED: Reducing capacity to `C = 3`, the resulting routes are:

```julia
Route 1: 1 => 5 => 1
Route 2: 1 => 4 => 1
Route 3: 1 => 3 => 1
Route 4: 1 => 2 => 1
Route 5: 1 => 1

```

---

<div class="post-metadata">

**Author:** ![Victor\_Neto](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/victor_neto/32/204390_2.png) [@Victor\_Neto](https://discourse.julialang.org/u/Victor_Neto)\
**Post date:** [November 29, 2023, 10:57pm UTC](https://discourse.julialang.org/t/subtour-elimination-constraints/106902/14 "2023-11-29T22:57:08Z")

</div>

Guys, I had taken a break because my head was no longer thinking lol, I’m coming back now, I’m going to do the tests soon.

---

<div class="post-metadata">

**Author:** ![Victor\_Neto](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/victor_neto/32/204390_2.png) [@Victor\_Neto](https://discourse.julialang.org/u/Victor_Neto)\
**Post date:** [November 30, 2023, 1:46am UTC](https://discourse.julialang.org/t/subtour-elimination-constraints/106902/15 "2023-11-30T01:46:39Z")

</div>

Hi guys, I did the tests with a toy model and got good results, the routes were well defined, I checked the values and they all matched the solution.  
I made a small change to these two restrictions to ensure that a bus does not pass on a route where there are no students.

# Flow in and out of other nodes is 1

[i in V; i != 1], sum(x[i, :, :]) \<= 1  
[i in V; i != 1], sum(x[:, i, :]) \<= 1

I would like to thank you immensely for your help.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [November 30, 2023, 1:58am UTC](https://discourse.julialang.org/t/subtour-elimination-constraints/106902/16 "2023-11-30T01:58:34Z")

</div>

Glad to hear. It can be tricky to get the MTZ and flow constraints right once you depart from a vanilla TSP.

The JuMP docs have an example of callback and iterative subtour elimination if the formulation begins to struggle:

> **[Traveling Salesperson Problem · JuMP](https://jump.dev/JuMP.jl/stable/tutorials/algorithms/tsp_lazy_constraints/)**
>
> Documentation for JuMP.

You’d likewise have to adapt it to look at all tours across `k`, but you should have a reasonable understanding now.