Subtour elimination constraints

Specific Domains Optimization (Mathematical)
, ,
1

Dear colleagues,
I am studying the school bus routing problem: Redirecting.
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?

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

2

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.

3

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.

4

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

5

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)

6

https://www.sciencedirect.com/science/article/abs/pii/S0377221713001586

7

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.

8

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.

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.

2 Likes
9

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.

10

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

V is a set of potentials stops

11

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

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
1 Like
12

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

13

Running odow’s code, I get:

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:

Route 1: 1 => 5 => 1
Route 2: 1 => 4 => 1
Route 3: 1 => 3 => 1
Route 4: 1 => 2 => 1
Route 5: 1 => 1
1 Like
14

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.

1 Like
15

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.

2 Likes
16

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:

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