Decompose twice at the same level.

[branchprediction]

Is your feature request related to a problem? Please describe.
According to the response by @guimarqu in #557 :

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

I attempted to do this with the following code which finds the correct optimal solution with SolveIpForm and no decomposition:

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.TreeSearchAlgorithm(explorestrategy = Coluna.Algorithm.BestDualBoundStrategy())
        # 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 axes
@axis(M_odd, 1:2:number_of_machines)
@axis(M_even, 2:2:number_of_machines)
# Define the other indices
M = 1: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)

# Assignment constraint is on all machines
@constraint(model, c2[j in J], sum(x[m, j] for m in 1:number_of_machines) == 1)

# Odd Machines
@constraint(model, c3[m in M_odd, 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, c5[m in M_odd, 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_odd, 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_odd, j in J, q in Qsub], b[m, j, q+1] - b[m, j, q] >= 0)
@constraint(model, c8[m in M_odd, j in J, q in Q], b[m, j, q] - x[m, j] <= 0)
@constraint(model, c9[m in M_odd, 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)
# Define makespan in master problem
@constraint(model, c4[m in 1:2: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))

# Even Machines
@constraint(model, c3_even[m in M_even, 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, c5_even[m in M_even, 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_even[m in M_even, 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_even[m in M_even, j in J, q in Qsub], b[m, j, q+1] - b[m, j, q] >= 0)
@constraint(model, c8_even[m in M_even, j in J, q in Q], b[m, j, q] - x[m, j] <= 0)
@constraint(model, c9_even[m in M_even, 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])
# Define makespan in master problem
@constraint(model, c4_even[m in 2:2: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))

@objective(model, Min, makespan)

@dantzig_wolfe_decomposition(model, decomposition_odd, M_odd)
@dantzig_wolfe_decomposition(model, decomposition_even, M_even)

subproblems_odd = getsubproblems(decomposition_odd)
subproblems_even = getsubproblems(decomposition_even)

specify!.(subproblems_odd, lower_multiplicity = 5, upper_multiplicity = 5)
specify!.(subproblems_even, lower_multiplicity = 5, upper_multiplicity = 5)

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

Which produces the following error:

ERROR: LoadError: Cannot decompose twice at the same level.
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] set_decomposition_tree!(::Model, ::Type{BlockDecomposition.DantzigWolfe}, ::BlockDecomposition.Axis{:M_even,Int64}) at /Path/To/BlockDecomposition/.../src/tree.jl:88
 [3] decompose_leaf(::Model, ::Type{BlockDecomposition.DantzigWolfe}, ::BlockDecomposition.Axis{:M_even,Int64}) at /Path/To/BlockDecomposition/.../src/tree.jl:134
 [4] top-level scope at /Path/To/BlockDecomposition/.../src/tree.jl:172
in expression starting at /Path/To/model.jl:82

Describe the feature or solution you'd like
This problem has 2 groups of identical sub-problems:

1. Machines with odd indices are identical to one another.
2. Machines with even indices are identical to one another.

This still falls within the block-diagonal structure needed to perform Dantzig-Wolfe decomposition, but it isn't clear whether it is currently possible to implement this within Coluna.

Do you feel able to code this feature or solution?
No.

[guimarqu]

You should create an axis with two entries, one for each type of subproblem :

@axis(M, 1:2)  # 1 for odd machines & 2 for even machines

You define your model as you did but each subproblem is defined only once.

Then, you declare the decomposition and specify the multiplicity of the two subproblems.

@dantzig_wolfe_decomposition(model, decomposition, M)

subproblems = getsubproblems(decomposition)

specify!(subproblems[1], lower_multiplicity = 5, upper_multiplicity = 5) # subproblem for odd machines
specify!(subproblems[2], lower_multiplicity = 5, upper_multiplicity = 5) # subproblem for even machines

Note that when you multiplicity, you cannot optimize the original formulation with a commercial solver anymore. The original formulation does not support the multiplicity, only the Dantzig-Wolfe reformulation supports it.

[branchprediction]

If I understand correctly, that would work for the case where sub-problem group 1 was identical to sub-problem group 2, but sub-problem group 1 (odd-indexed machines) has an additional constraint (c15). The two requirements are:

1. The axis range S must occur in the index definition of each constraint.
2. Some constraints only belong to a subset of the indices of the axis range.

Example: In my scheduling problem, constraint c15 should only be generated when s=1, but S=1:2 and therefore I need to somehow "do nothing" when s=2.

I made an attempt to accommodate the two requirements by making the machine index m conditional upon the value of s.
It works correctly except for c15 where @dantzig_wolfe_decomposition appears to dislike the fact that nothing happens when s=1. Removing c15 works for both "SolveIpForm and no decomposition" as well as with Dantzig-Wolfe Decomposition.

Any ideas how I might accommodate sub-problem group 1 having the additional c15 constraint?

UPDATE: Actually, I think this seems to run properly for both "SolveIpForm and no decomposition" as well as with Dantzig-Wolfe Decomposition. At the very least it runs the column generation, even if it isn't producing a result.

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_subproblem_groups = 2
number_of_machines = 10
number_of_jobs = 21
number_of_nonavailability_periods_per_machine = 1
# Define the decomposition axis
@axis(S, 1:number_of_subproblem_groups)  # s=1 indexes odd machines, s=2 indexes even machines, s in S.
# Define the other indices
M = 1:number_of_machines  # indices of all machines
M_odd = 1:2:number_of_machines  # indices of odd machines
M_even = 2:2:number_of_machines  # indices of even machines
# In constraints, we need machine indices conditioned upon the value of s in S, therefore need a ragged array.
MV = Array{Array{Int}}(undef, 2)  # Can be a ragged array, therefore good for differing numbers of sub-problems per sub-problem group.
MV[1] = M_odd  # MV[s] gives odd_machine_indices when s==1
MV[2] = M_even  # MV[s] gives even_machine_indices when s==2
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 = 1000.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 1:number_of_machines) == 1)
@constraint(model, c3[s in S, m in MV[s], 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])
# IMPORTANT: Using the axis variable in the second argument of constraint results in it being put in the sub-problem.
# Use the explicit definition of the axis variable if you want it to be put in the master problem.
for q in Q) - makespan <= - sum(Q_durs[m, q] for q in Q))
@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))
@constraint(model, c5[s in S, m in MV[s], 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[s in S, m in MV[s], 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[s in S, m in MV[s], j in J, q in Qsub], b[m, j, q+1] - b[m, j, q] >= 0)
@constraint(model, c8[s in S, m in MV[s], j in J, q in Q], b[m, j, q] - x[m, j] <= 0)
@constraint(model, c9[s in S, m in MV[s], 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])

# All the above constraints (except c2 and c4 which are in the master problem) occurred in both sub-problem groups.
# However, c15 only occurs in sub-problem group 1 (machines with odd indices).  Therefore, we need a ragged array 
# with no indices in the first element.  In the most general case, we would need a ragged array for each @constraint macro.
MVc15 = Array{Array{Int}}(undef, 2)
MVc15[1] = M_odd
MVc15[2] = []
@constraint(model, c15[s in S, m in MVc15[s]], sum(x[m, j] for j in J) + y[m, 1] >= 1)

@objective(model, Min, makespan)

@dantzig_wolfe_decomposition(model, decomposition, S)

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

specify!(subproblems[1], lower_multiplicity = 0, upper_multiplicity = 5)
specify!(subproblems[2], lower_multiplicity = 0, upper_multiplicity = 5)

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

[guimarqu]

Yes but you don't need to duplicate the constraints for each machine of each type otherwise you don't take advantage of the multiplicity (here, you still have 5 subproblems for each type of machines but you solve them at the same time).

The goal of multiplicity is to have only one subproblem for each type of machine. A subproblem will generate assignement pattern for a type of machine. Then the multiplicity allows the master to select from 0 to 5 assignment patterns generated by the subproblem.

You can remove MV[s] in all constraints and if c15 is only in the second subproblem, you can write :

@constraint(model, c15[S[2]], sum(x[s, j] for j in J) + y[s, 1] >= 1)

[branchprediction]

Ok, I think I understand now. I'm guessing that the variables still get indexed by machine and that the change of index from machines to sub-problem group index only applies to the sub-problem constraints, and therefore master constraints still remain indexed by machine. Please advise whether the variables should also be indexed by s instead of m and how that would change the master constraints. Right now it seems to run but I can no longer use SolveIpForm to check the solution; which makes sense and matches what you indicated.

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()  # No longer an option for this multiplicity formulation
    ),
    "default_optimizer" => Gurobi.Optimizer # Gurobi for the master & the sub-problems
)

number_of_subproblem_groups = 2
number_of_machines = 10
number_of_jobs = 21
number_of_nonavailability_periods_per_machine = 1
# Define the decomposition axis
@axis(S, 1:number_of_subproblem_groups)  # s=1 indexes odd machines, s=2 indexes even machines, s in S.
# Define the other indices
M = 1:number_of_machines  # indices of all 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 = 1000.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 1:number_of_machines) == 1)  # Master constraint, therefore still indexed by machine index.

# @constraint(model, c3[s in S, m in MV[s], 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, c3[s in S, q in Q], sum(J_durs[j] * x[s, j] for j in J) + sum(alpha * w[s, k] for k in 1:(q-1)) + big_m * y[s, q] <= - sum(Q_durs[s, k] for k in 1:(q-1)) + big_m + Q_starts[s, q])

@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))  # Master constraint, therefore still indexed by machine index.

# @constraint(model, c5[s in S, m in MV[s], 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, c5[s in S, q in Q], Q_starts[s, q] * y[s, q] + sum(J_durs[j] * b[s, j, q] for j in J) + sum(alpha * w[s, k] for k in 1:(q-1)) + w[s, q] >= Q_starts[s, q] - sum(Q_durs[s, k] for k in 1:(q-1)))

# @constraint(model, c6[s in S, m in MV[s], 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, c6[s in S, q in Q], sum(J_durs[j] * b[s, j, q] for j in J) + sum(alpha * w[s, k] for k in 1:(q-1)) + w[s, q] <= Q_starts[s, q] - sum(Q_durs[s, k] for k in 1:(q-1)))

# @constraint(model, c7[s in S, m in MV[s], j in J, q in Qsub], b[m, j, q+1] - b[m, j, q] >= 0)
@constraint(model, c7[s in S, j in J, q in Qsub], b[s, j, q+1] - b[s, j, q] >= 0)

# @constraint(model, c8[s in S, m in MV[s], j in J, q in Q], b[m, j, q] - x[m, j] <= 0)
@constraint(model, c8[s in S, j in J, q in Q], b[s, j, q] - x[s, j] <= 0)

# @constraint(model, c9[s in S, m in MV[s], 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, c9[s in S, q in Qsub], sum(J_durs[j] * b[s, j, q+1] for j in J) - sum(J_durs[j] * b[s, j, q] for j in J) - (1.0 - alpha) * w[s, q] <= Q_starts[s, q+1] - Q_ends[s, q])


# All the above constraints (except c2 and c4 which are in the master problem) occurred in both sub-problem groups.
# However, c15 only occurs in sub-problem group 1 (machines with odd indices).  Therefore, the API allows us to
# subset the axis variable S[1].
@constraint(model, c15[s in S[1]], sum(x[s, j] for j in J) + y[s, 1] >= 1)

@objective(model, Min, makespan)

@dantzig_wolfe_decomposition(model, decomposition, S)

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

specify!(subproblems[1], lower_multiplicity = 0, upper_multiplicity = 5)
specify!(subproblems[2], lower_multiplicity = 0, upper_multiplicity = 5)

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

[guimarqu]

Ok, I think I understand now. I'm guessing that the variables still get indexed by machine and that the change of index from machines to sub-problem group index only applies to the sub-problem constraints, and therefore master constraints still remain indexed by machine.

No variables and constraints, even in the master, should be indexed by the machine types.

@variable(model, x[s in S, j in J], Bin)  # equals one if job j is assigned to machine m
...

Please advise whether the variables should also be indexed by s instead of m and how that would change the master constraints. Right now it seems to run but I can no longer use SolveIpForm to check the solution; which makes sense and matches what you indicated.

Yes documentation must be improved, and some fixes must be done for the case with identical subproblems.

[branchprediction]

Ok, then I suppose the last question is how would the following master constraints change and still give the same meaning?

@constraint(model, c2[j in J], sum(x[m, j] for m in 1:number_of_machines) == 1)
@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))

If x is no longer indexed by m, but rather s, am I correct that they would become the following?

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

I guess the way I would have to think about "s in S" is not as "s in 1:2" but rather "s in (Columns[1] union Columns[2])" where Columns[sp] are the columns generated by identical sub-problem group sp.

[guimarqu]

Yes, that's it !