Вычисление несводимой несогласованной подсистемы (IIS) с использованием Julia JuMP (Gurobi)

Пытаюсь вычислить IIS для моей до глупости чрезмерно сложной модели.

Я включу всю модель для ясности:

using JuMP
using Gurobi
import XLSX
roster = Model(Gurobi.Optimizer)

Intern = 1:11 #i
Week = 1:52  #k
Rotation = 1:23 #j
Leave_week = 1:3
Dec_leave = 1:2
M = 1000
clins = 7:52
non_clins = 5:52
early = 5:28
gen = [5,8]

@variables(roster, begin
    x[Intern,Week, Rotation], Bin
    y[Intern,Week, Rotation], Bin
    L[Leave_week, Week], Bin
    D[Dec_leave, Intern], Bin
    s[Intern, Week], Bin
    g[Intern, gen], Bin
end
)

#physical constraint
@constraint(roster, phys[i in Intern, k in Week],
    sum(x[i,k,j] for j in Rotation) == 1)

#Rotation capacity
rots    = [1,2,3,4,5,6,7,8,9,10,12,13,15,16,17,18,19]
cap_rhs = [2,1,1,1,1,1,1,1,1, 1, 1, 1, 2, 1, 1, 2, 1, 1]
cap = @constraint(roster, [(b,d) in zip(rots, cap_rhs), k in 1:52],
    sum(x[i,k,b] for i in Intern) <= d)

#dispensary
disp  = @constraint(roster, [i in Intern], sum(x[i,k,j] for k in Week, j in 14:18) >= 5)
disp1 = @constraint(roster, [i in Intern], sum(x[i,k,j] for k in 29:40, j in 14:18) >= 1)
disp2 = @constraint(roster, [i in Intern], sum(x[i,k,j] for k in 41:52, j in 14:18) >= 1)
clay_cap_o = @constraint(roster, [k in 1:4],
        sum(x[i,k,14] for i in Intern) <=3)
clay_cap_o = @constraint(roster, [k in 5:52],
        sum(x[i,k,14] for i in Intern) <=2)

#Orientation
IP_1 = @constraint(roster, [i in Intern], sum(x[i,k,1] for k in 1:6) >= 1)
IP
The C API of Gurobi.jl has been rewritten to expose the complete C API, and
all old functions have been removed. For more information, see the Discourse
announcement: https://discourse.julialang.org/t/ann-upcoming-breaking-changes-to-cplex-jl-and-gurobi-jl

Here is a brief summary of the changes.

...

1 = @constraint(roster, [i in Intern], sum(s[i,k] for k in 1:5) <= 1)
IP
The C API of Gurobi.jl has been rewritten to expose the complete C API, and
all old functions have been removed. For more information, see the Discourse
announcement: https://discourse.julialang.org/t/ann-upcoming-breaking-changes-to-cplex-jl-and-gurobi-jl

Here is a brief summary of the changes.

...

11 = @constraint(roster, sum(s[i,k] for i in Intern, k in 1:5) == 10)
IP
The C API of Gurobi.jl has been rewritten to expose the complete C API, and
all old functions have been removed. For more information, see the Discourse
announcement: https://discourse.julialang.org/t/ann-upcoming-breaking-changes-to-cplex-jl-and-gurobi-jl

Here is a brief summary of the changes.

...

2 = @constraint(roster, [i in Intern, k in 1:4], x[i,k,1] == s[i,k])
# IP_lazy = @constraint(roster, [(i,k) in zip(Intern, [1 1 2 2 3 3 4 4 5 5 6])], x[i,k,1] ==1)
orien = @constraint(roster, [i in Intern], sum(x[i,k,j] for k in 1:4, j in [1,14,15,16,17,18] ) == 4)
orien1 = @constraint(roster, [i in Intern, j in [1,14,15,16,17,18]], sum(x[i,k,j] for k in 1:4 ) <= 2)

#leave
# 2 weeks leave
@constraint(roster, [i in Intern],
                    sum(x[i,k,j] for k in Week, j in 20:22) == 2)
week1_dvar = @constraint(roster, sum(L[1, k] for k in 17:22) == 1)
@constraint(roster, week1[k in 17:22], sum(x[i,k,20] for i in Intern) == 11*L[1,k])
@constraint(roster, week2
The C API of Gurobi.jl has been rewritten to expose the complete C API, and
all old functions have been removed. For more information, see the Discourse
announcement: https://discourse.julialang.org/t/ann-upcoming-breaking-changes-to-cplex-jl-and-gurobi-jl

Here is a brief summary of the changes.

...

dvar[l in 2:3], sum(L[l, k] for k in 35:41) == 1)
@constraint(roster, week2_3[(l, j, rhs) in zip(2:3, 21:22, [6,5]), k in 35:41],
    sum(x[i,k,j] for i in Intern) == rhs*L[l, k] )
@constraint(roster, max_leave[i in Intern], sum(x[i,k,j] for j in 20:22, k in Week) ==2)
## - Dec_leave
@constraint(roster,[i in Intern], sum(D[l,i] for l in 1:2) == 1)
@constraint(roster, [(l,d) in zip(1:2,[6,5])], sum(D[l,i] for i in Intern) == d)
@constraint(roster, [i in Intern, (l,b) in zip(1:2, [49:50, 51:52])],
    sum(x[i,k,23] for k in b) == 2*D[l,i])
@constraint(roster, [i in Intern], sum(x[i,k,23] for k in Week) == 2)

#MIC
MIC_1_dvar = @constraint(roster, [i in Intern], sum(y[i,k,4] for k in 5:27 ) == 1)
MIC Gurobi.computeIIS(roster)dvar = @constraint(roster, [i in Intern], sum(y[i,k,4] for k in 29:51 ) == 1)
MIC = @constraint(roster, [ i in Intern, k in 5:27],
    2 - sum(x[i, k + alpha, 4] for alpha in 0:1 ) <= M*(1-y[i,k,4]))
MIC = @constraint(roster, [ i in Intern, k in 29:51],
    2 - sum(x[i, k + alpha, 4] for alpha in 0:1 ) <= M*(1-y[i,k,4]))

#gen_med
g_vars = @constraint(roster, [i in Intern], sum(g[i,m] for m in gen) ==1)
gen_duration_dvar = @constraint(roster, [(b,d) in zip(gen,[6,7]), i in Intern],
                sum(y[i,k,b] for k in 1:(52 - (d-1) ) ) == g[i,b])
gen_limit = @constraint(roster, [(b,d) in zip(gen,[6,7]), i in Intern],
                sum(x[i,k,b] for k in Week) == g[i,b]*d)
gen_durations = @constraint(roster, [(b,d) in zip(gen,[6,7]),
                i in Intern, k in 1:(52 - (d-1) )],
                d - sum(x[i, k + alpha, b] for alpha in 0:(d-1) ) <= M*(1-y[i,k,b]))
ed_with_gen = @constraint(roster, [i in Intern, k in 2:50], y[i,k,23] - x[i,k-1,8] - x[i,k+2,8] <= (1-g[i,5]))

#qum
qum_1 = @constraint(roster, [i in Intern], sum(x[i,k,13] for k in early) >= 1)
qum_2 = @constraint(roster, [i in Intern], sum(x[i,k,13] for k in 1:39) == 2)

# duration
dur_rot =   [2,6,7,9,10,11,19]
durs =      [2,4,2,3, 3, 4, 2]
duration_dvar = @constraint(roster, [(b,d) in zip(dur_rot, durs), i in Intern],
                sum(y[i,k,b] for k in 1:(52 - (d-1) ) ) == 1)
durations = @constraint(roster, [(b,d) in zip(dur_rot, durs),
                i in Intern, k in 1:(52 - (d-1) )],
                d - sum(x[i, k + alpha, b] for alpha in 0:(d-1) ) <= M*(1-y[i,k,b]))
                
AP_dur_var = @constraint(roster, [i in Intern], sum(y[i,k,3] for k in 5:35) == 1)
AP_dur     = @constraint(roster, [i in Intern, k in 5:35],
                2 - sum(x[i,k + alpha, 3] for alpha in 0:1) <= M*(1 - y[i,k,3]))
AP_third   = @constraint(roster,[i in Intern], sum(x[i,k,3] for k in 37:52) == 1)



# rotations_lengths
completion = @constraint(roster,
    [(j,c,d) in zip([1,2,3,4,6,7,9,10,11,12,19],
    [ 1:28,clins,non_clins, non_clins,clins,clins, clins, clins, 29:52, early, clins],
                    [3,2,3,4,4,2,3, 3, 4, 1, 2]), i in Intern],
                  sum(x[i,k,j] for k in c) == d)
IP_soft = @constraint(roster, [i in Intern], sum(x[i,k,1] for k in Week) >= 5)
whole_year = @constraint(roster, [(j,d) in zip([1,2,3,4,6,7,9,10,11,12,13,19],
                                                [5,2,3,4,4,2,3, 3, 4, 1, 2, 2]), i in Intern],
                                                sum(x[i,k,j] for k in Week) == d)

#  public holiday constraints
no_pubs = @constraint(roster,
    [i in Intern, k in [4,7,10,13,14,22,24,29,39,44,52], j in [12,13]],
    x[i,k,j] == 0 )

z = @expression(roster, sum(x[i,k,j] for i in Intern, j in Rotation, k in Week))
obj_z = @objective(roster, Max, z)

optimize!(roster)

Пробовал несколько подходов, например рекомендованные здесь и здесь, но возникли ошибки.

Пытался Gurobi.computeIIS(roster), но получилось:

The C API of Gurobi.jl has been rewritten to expose the complete C API, and
all old functions have been removed. For more information, see the Discourse
announcement: https://discourse.julialang.org/t/ann-upcoming-breaking-changes-to-cplex-jl-and-gurobi-jl

Here is a brief summary of the changes.

...

Любая помощь или совет приветствуются. Пожалуйста, не стесняйтесь редактировать пример для ясности.

Большое спасибо.


julia gurobi julia-jump integer-programming
person Chris Swan    schedule 12.11.2020    source источник
comment
Я знаю, что мои переменные Решения ошибочны - планирую обновить их в будущем, чтобы избежать всех ужасных дваров Big M. Но, на мой взгляд, это должно просто уничтожить время выполнения, а не обязательно сделать его невозможным.   -  person Chris Swan    schedule 12.11.2020

Ответы (1)


arrow_upward
1
arrow_downward

Это немного продвинуто, и часть сантехники отсутствует (в настоящее время вы должны использовать direct_model), но вы можете пойти:

using JuMP, Gurobi
model = direct_model(Gurobi.Optimizer())
@variable(model, x >= 0)
@constraint(model, c1,     x <= -1)
@constraint(model, c2, 2 * x <=  1)
optimize!(model)
@assert termination_status(model) == MOI.INFEASIBLE
compute_conflict!(model)

julia> MOI.get(model, MOI.ConstraintConflictStatus(), LowerBoundRef(x))
IN_CONFLICT::ConflictParticipationStatusCode = 1

julia> MOI.get(model, MOI.ConstraintConflictStatus(), c1)
IN_CONFLICT::ConflictParticipationStatusCode = 1

julia> MOI.get(model, MOI.ConstraintConflictStatus(), c2)
NOT_IN_CONFLICT::ConflictParticipationStatusCode = 0
person Oscar Dowson    schedule 12.11.2020
comment
На самом деле, когда я оптимизировал с помощью direct_model, он нашел оптимальное решение!?! - person Chris Swan; 13.11.2020