Error in TreeSearchAlgorithm but not SolveIpForm for MILP

[branchprediction]

Describe the bug
Simple decomposable MILP scheduling problem solves with solver = Coluna.Algorithm.SolveIpForm() in 1 second but fails with solver = Coluna.Algorithm.TreeSearchAlgorithm().

To Reproduce

using JuMP, BlockDecomposition, Coluna, Gurobi

optimizer = optimizer_with_attributes(
    Coluna.Optimizer,
    "params" => Coluna.Params(
        # solver = Coluna.Algorithm.ColumnGeneration()
        solver = Coluna.Algorithm.TreeSearchAlgorithm() # default BCP
        # solver = Coluna.Algorithm.SolveIpForm()
    ),
    "default_optimizer" => Gurobi.Optimizer # Gurobi for the master & the sub-problems
)

number_of_machines = 10
number_of_jobs = 21
number_of_nonavailability_periods_per_machine = 1
# Define the decomposition axis
@axis(M, 1:number_of_machines)
# Define the other indices
M_odd = 1:2:number_of_machines
J = 1:number_of_jobs
Q = 1:number_of_nonavailability_periods_per_machine
Qsub = 1:(number_of_nonavailability_periods_per_machine-1)
# Define parameters
big_m = 100000.0
alpha = 0.5
Q_starts = [0.0; 20.0; 0.0; 20.0; 0.0; 20.0; 0.0; 20.0; 0.0; 20.0]
Q_ends = [2.0; 22.0; 2.0; 22.0; 2.0; 22.0; 2.0; 22.0; 2.0; 22.0]
Q_durs = [2.0; 2.0; 2.0; 2.0; 2.0; 2.0; 2.0; 2.0; 2.0; 2.0]
J_durs = [19.0, 19.0, 18.0, 18.0, 17.0, 17.0, 16.0, 16.0, 15.0, 15.0, 14.0, 14.0, 13.0, 13.0, 12.0, 12.0, 11.0, 11.0, 10.0, 10.0, 10.0]

model = BlockModel(optimizer)
@variable(model, x[m in M, j in J], Bin)
@variable(model, y[m in M, q in Q], Bin)
@variable(model, b[m in M, j in J, q in Q], Bin)
@variable(model, 0.0<=w[m in M, q in Q]<=1000.0)
@variable(model, 0.0<=makespan<=1000.0)

@constraint(model, c2[j in J], sum(x[m, j] for m in M) == 1)
@constraint(model, c3[m in M, q in Q], sum(J_durs[j] * x[m, j] for j in J) + sum(alpha * w[m, k] for k in 1:(q-1)) + big_m * y[m, q] <= - sum(Q_durs[m, k] for k in 1:(q-1)) + big_m + Q_starts[m, q])
@constraint(model, c4[m in M], sum(J_durs[j] * x[m, j] for j in J) - sum(Q_durs[m, q] * y[m, q] for q in Q) + sum(alpha * w[m, q] for q in Q) - makespan <= - sum(Q_durs[m, q] for q in Q))
@constraint(model, c5[m in M, q in Q], Q_starts[m, q] * y[m, q] + sum(J_durs[j] * b[m, j, q] for j in J) + sum(alpha * w[m, k] for k in 1:(q-1)) + w[m, q] >= Q_starts[m, q] - sum(Q_durs[m, k] for k in 1:(q-1)))
@constraint(model, c6[m in M, q in Q], sum(J_durs[j] * b[m, j, q] for j in J) + sum(alpha * w[m, k] for k in 1:(q-1)) + w[m, q] <= Q_starts[m, q] - sum(Q_durs[m, k] for k in 1:(q-1)))
@constraint(model, c7[m in M, j in J, q in Qsub], b[m, j, q+1] - b[m, j, q] >= 0)
@constraint(model, c8[m in M, j in J, q in Q], b[m, j, q] - x[m, j] <= 0)
@constraint(model, c9[m in M, q in Qsub], sum(J_durs[j] * b[m, j, q+1] for j in J) - sum(J_durs[j] * b[m, j, q] for j in J) - (1.0 - alpha) * w[m, q] <= Q_starts[m, q+1] - Q_ends[m, q])
@constraint(model, c15[m in M_odd], sum(x[m, j] for j in J) + y[m, 1] >= 1)

@objective(model, Min, makespan)

@dantzig_wolfe_decomposition(model, decomposition, M)

# master = getmaster(decomposition)
# subproblems = getsubproblems(decomposition)

# specify!.(subproblems, lower_multiplicity = 1, upper_multiplicity = 1)

optimize!(model)
for m in M
    a = collect(value.(x[m, j]) for j in J)
    println(a)
end

Expected behavior
This problem has a known optimal solution of 33. Expected both SolveIpForm and TreeSearchAlgorithm to produce the same result. Instead got the following:

Result using solver = Coluna.Algorithm.SolveIpForm() :

Coluna
Version 0.3.9 | 2021-05-15 | https://github.com/atoptima/Coluna.jl
┌ Warning: No initial primal bound and no cost for global artificial variables.
│ Cost of global artificial variables set to 100000.0
└ @ Coluna /Path/To/optimize.jl:20
┌ Warning: No initial primal bound and no cost for local artificial variables.
│ Cost of local artificial variables set to 10000.0
└ @ Coluna /Path/To/optimize.jl:33
Found primal solution of 33.0000 
 ──────────────────────────────────────────────────────────────────
                           Time                   Allocations      
                   ──────────────────────   ───────────────────────
 Tot / % measured:      1197s / 0.03%            355MiB / 0.38%    

 Section   ncalls     time   %tot     avg     alloc   %tot      avg
 ──────────────────────────────────────────────────────────────────
 Coluna         1    351ms   100%   351ms   1.35MiB  100%   1.35MiB
 ──────────────────────────────────────────────────────────────────
[ Info: Terminated
[ Info: Primal bound: 33.0
[ Info: Dual bound: 33.0
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0]
[0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
[0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

Result using solver = Coluna.Algorithm.TreeSearchAlgorithm():

Coluna
Version 0.3.9 | 2021-05-15 | https://github.com/atoptima/Coluna.jl
┌ Warning: No initial primal bound and no cost for global artificial variables.
│ Cost of global artificial variables set to 100000.0
└ @ Coluna /Path/To/optimize.jl:20
┌ Warning: No initial primal bound and no cost for local artificial variables.
│ Cost of local artificial variables set to 10000.0
└ @ Coluna /Path/To/optimize.jl:33
***************************************************************************************
**** BaB tree root node
**** Local DB = -Inf, global bounds : [ -Inf , Inf ], time = 0.00 sec.
***************************************************************************************
  <it=  1> <et= 0.01> <mst= 0.00> <sp= 0.01> <cols=10> <al= 0.00> <DB=60000.0000> <mlp=200000.0000> <PB=Inf>
  <it=  2> <et= 0.02> <mst= 0.00> <sp= 0.01> <cols= 0> <al= 0.00> <DB=100000.0000> <mlp=100000.0000> <PB=Inf>
[ Info: Column generation algorithm has converged.
# <it=  1> <et= 0.03> <mst= 0.00> <sp= 0.01> <cols= 0> <al= 0.00> <DB=100000.0000> <mlp=100000.0000> <PB=Inf>
[ Info: Phase one determines infeasibility.
[ Info: No branching candidates found. No children will be generated.
 ────────────────────────────────────────────────────────────────────────────────────
                                             Time                   Allocations      
                                     ──────────────────────   ───────────────────────
          Tot / % measured:               24.4s / 0.11%           84.4MiB / 4.77%    

 Section                     ncalls     time   %tot     avg     alloc   %tot      avg
 ────────────────────────────────────────────────────────────────────────────────────
 Coluna                           1   27.6ms   100%  27.6ms   4.02MiB  100%   4.02MiB
   SolveLpForm                    3   2.89ms  10.5%   965μs    227KiB  5.51%  75.6KiB
   Update reduced costs           3    473μs  1.71%   158μs    655KiB  15.9%   218KiB
   Inserting columns              3    162μs  0.59%  54.1μs   66.5KiB  1.61%  22.2KiB
   Cleanup columns                3   36.2μs  0.13%  12.1μs   1.36KiB  0.03%     464B
   Restore/remove records         2   25.3μs  0.09%  12.6μs      352B  0.01%     176B
   Getting primal solution        3   20.1μs  0.07%  6.71μs      272B  0.01%    90.7B
   Store records                  1   17.4μs  0.06%  17.4μs      592B  0.01%     592B
   Update Lagrangian bound        3   6.35μs  0.02%  2.12μs     0.00B  0.00%    0.00B
   Smoothing update               3   3.29μs  0.01%  1.10μs     48.0B  0.00%    16.0B
 ────────────────────────────────────────────────────────────────────────────────────
[ Info: Terminated
[ Info: Primal bound: Inf
[ Info: Dual bound: Inf
┌ Warning: Coluna did not find a primal feasible solution.
└ @ Coluna /Path/To/MOIwrapper.jl:651
┌ Warning: Coluna did not find a primal feasible solution.
└ @ Coluna /Path/To/MOIwrapper.jl:651
┌ Warning: Coluna did not find a primal feasible solution.
└ @ Coluna /Path/To/MOIwrapper.jl:651
.
.
.
└ @ Coluna /Path/To/MOIwrapper.jl:651
┌ Warning: Coluna did not find a primal feasible solution.
└ @ Coluna /Path/To/MOIwrapper.jl:651
┌ Warning: Coluna did not find a primal feasible solution.
└ @ Coluna /Path/To/MOIwrapper.jl:651
┌ Warning: Coluna did not find a primal feasible solution.
└ @ Coluna /Path/To/MOIwrapper.jl:651
[NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN]
┌ Warning: Coluna did not find a primal feasible solution.
└ @ Coluna /Path/To/MOIwrapper.jl:651
.
.
.
etc

Environment (please complete the following information):

• Julia Version 1.4.1
• OS: Linux (x86_64-linux-gnu)
• BlockDecomposition v1.3.0
• Coluna v0.3.9
• Gurobi v0.9.14
• JuMP v0.21.8

Additional context
This problem has a known optimal solution of 33. Expected dantzig_wolfe_decomposition to result in a master problem containing constraints c2 (assignment of jobs to machines) and c4 (defines makespan >= machine makespans) and other constraints in the sub-problems. The makespan variable is master-only, and it's the only variable with a non-zero objective coefficient. All machines are independent similar to GAP example, so I was really expecting this to work. Furthermore, there are 2 sets of identical sub-problems: machines with odd indices (there are 5 of them) are identical to one another, and machines with even indices (there are 5 of them) are identical to one another. Not sure whether I need to duplicate constraints for M_odd, M_even instead of just M and use lower_multiplicity and upper_multiplicity to make this work. If so, what values to use for multiplicities (5, 5?). Also, that would require multiple @axis and presumably multiple @dantzig_wolfe_decomposition and I'm not sure if that is supported. Do I have to mark the sub-problems as identical some how or was that just in earlier versions of Coluna? Please advise.

[guimarqu]

Hi,

Constraint c4 is assign to subproblems because you use the axis M in its definition. You should do :

@constraint(model, c4[m in 1:number_of_machines], sum(J_durs[j] * x[m, j] for j in J) - sum(Q_durs[m, q] * y[m, q] for q in Q) + sum(alpha * w[m, q] for q in Q) - makespan <= - sum(Q_durs[m, q] for q in Q))

and the constraint will go into the master.

You can check that with method BlockDecomposition.annotation(model, c4[1]) :

julia> BlockDecomposition.annotation(model, c4[1])
Annotation(BlockDecomposition.Master, BlockDecomposition.DantzigWolfe, lm = 1.0, um = 1.0, id = 2)

Furthermore, there are 2 sets of identical sub-problems: machines with odd indices (there are 5 of them) are identical to one another, and machines with even indices (there are 5 of them) are identical to one another. Not sure whether I need to duplicate constraints for M_odd, M_even instead of just M and use lower_multiplicity and upper_multiplicity to make this work. If so, what values to use for multiplicities (5, 5?). Also, that would require multiple @axis and presumably multiple @dantzig_wolfe_decomposition and I'm not sure if that is supported. Do I have to mark the sub-problems as identical some how or was that just in earlier versions of Coluna? Please advise.

You can use lower multiplicity and upper multiplicity. (5,5) if the solution should use the 5 machines or (0,5) if some machines can stay unused. Both are equivalent if the subproblem can generate a solution that does not assign any job to a machine. However, multiplicity is not fully supported at the moment :

The current issue that in some cases a solution may be declared integer (because all original variables have integer values), but no disaggregated integer solution exist. This should be addressed in the future. As I understand, this may happen only if there are general integer subproblem variables and subproblems multiplicity is greater than one (although I am not 100% sure).

Originally posted by @rrsadykov in #551 (comment)

This is also a problem for branching because when all original variables are integer, we cannot branch anymore. The algorithm then stops prematurely.

[guimarqu]

Coluna should throw an error because the master variable makespan should not be in the subproblems.

[guimarqu]

Hi,

With the latest version of BlockDecomposition, the following error is now thrown:

ERROR: LoadError: BlockDecomposition.MasterVarInDwSp(makespan, c4[1] : 19 x[1,1] + 19 x[1,2] + 18 x[1,3] + 18 x[1,4] + 17 x[1,5] + 17 x[1,6] + 16 x[1,7] + 16 x[1,8] + 15 x[1,9] + 15 x[1,10] + 14 x[1,11] + 14 x[1,12] + 13 x[1,13] + 13 x[1,14] + 12 x[1,15] + 12 x[1,16] + 11 x[1,17] + 11 x[1,18] + 10 x[1,19] + 10 x[1,20] + 10 x[1,21] - 2 y[1,1] + 0.5 w[1,1] - makespan ≤ -2.0)

Subproblem multiplicity is indicated as an alpha feature in the documentation. Progress will be made soon.

Thanks for the report.