Identical subproblems
Let us see an example of resolution using the advantage of identical subproblems with Dantzig-Wolfe and a variant of the Generalized Assignment Problem.
Consider a set of machine type T = 1:nb_machine_types and a set of jobs J = 1:nb_jobs. A machine type t has a resource capacity Q[t] and the factory contains U[t] machines of type t. A job j assigned to a machine of type t has a cost c[t,j] and consumes w[t,j] resource units of the machine of type t.
Consider the following instance :
nb_machine_types = 2;
nb_jobs = 8;
J = 1:nb_jobs;
Q = [10, 15];
U = [3, 2]; # 3 machines of type 1 & 2 machines of type 2
c = [10 11 13 11 12 14 15 8; 20 21 23 21 22 24 25 18];
w = [4 4 5 4 4 3 4 5; 5 5 6 5 5 4 5 6];
#Here is the JuMP model to optimize this instance with a classic solver :
using JuMP, GLPK;
T1 = [1, 2, 3]; # U[1] machines
T2 = [4, 5]; # U[2] machines
M = union(T1, T2);
m2t = [1, 1, 1, 2, 2]; # machine id -> type id
model = Model(GLPK.Optimizer);
@variable(model, x[M, J], Bin); # 1 if job j assigned to machine m
@constraint(model, cov[j in J], sum(x[m,j] for m in M) == 1);
@constraint(model, knp[m in M], sum(w[m2t[m],j] * x[m,j] for j in J) <= Q[m2t[m]]);
@objective(model, Min, sum(c[m2t[m],j] * x[m,j] for m in M, j in J));
optimize!(model);
objective_value(model)114.0You can decompose over the machines by defining an axis on M. However, if you want to take advantage of the identical subproblems, you must define the formulation as follows :
using BlockDecomposition, Coluna, JuMP, GLPK;
const BD = BlockDecomposition
coluna = optimizer_with_attributes(
Coluna.Optimizer,
"params" => Coluna.Params(
solver = Coluna.Algorithm.TreeSearchAlgorithm() # default BCP
),
"default_optimizer" => GLPK.Optimizer # GLPK for the master & the subproblems
);
@axis(T, 1:nb_machine_types);
model = BlockModel(coluna);
@variable(model, x[T, J], Bin); # 1 if job j assigned to machine m
@constraint(model, cov[j in J], sum(x[t,j] for t in T) == 1);
@constraint(model, knp[t in T], sum(w[t] * x[t,j] for j in J) <= Q[t]);
@objective(model, Min, sum(c[t,j] * x[t,j] for t in T, j in J));We assign jobs to a type of machine and we define one knapsack constraint for each type. This formulation cannot be solved as it stands with a commercial solver.
Then, we decompose and specify the multiplicity of each knapsack subproblem :
@dantzig_wolfe_decomposition(model, dec_on_types, T);
sps = getsubproblems(dec_on_types)
for t in T
specify!(sps[t], lower_multiplicity = 0, upper_multiplicity = U[t]);
end
getsubproblems(dec_on_types)2-element Vector{BlockDecomposition.SubproblemForm}:
Subproblem formulation for T = 1 contains : 0.0 <= multiplicity <= 3.0
Subproblem formulation for T = 2 contains : 0.0 <= multiplicity <= 2.0
We see that subproblem for machine type 1 has an upper multiplicity equals to 3, and the second subproblem for machine type 2 has an upper multiplicity equals to 2. It means that we can use at most 3 machines of type 1 and at most 2 machines of type 2.
We can then optimize
optimize!(model);Coluna
Version 0.8.1 | https://github.com/atoptima/Coluna.jl
***************************************************************************************
**** B&B tree root node
**** Local DB = -Inf, global bounds: [ -Inf , Inf ], time = 0.00 sec.
***************************************************************************************
<st= 1> <it= 1> <et= 0.00> <mst= 0.00> <sp= 0.00> <cols= 2> <al= 0.00> <DB=-39828.0000> <mlp=80000.0000> <PB=Inf>
<st= 1> <it= 2> <et= 0.00> <mst= 0.00> <sp= 0.00> <cols= 2> <al= 0.00> <DB=-69747.0000> <mlp=50059.0000> <PB=Inf>
<st= 1> <it= 3> <et= 0.00> <mst= 0.00> <sp= 0.00> <cols= 2> <al= 0.00> <DB=-89715.0000> <mlp=30082.0000> <PB=Inf>
<st= 1> <it= 4> <et= 0.00> <mst= 0.00> <sp= 0.00> <cols= 2> <al= 0.00> <DB=-79810.0000> <mlp=20131.0000> <PB=Inf>
<st= 1> <it= 5> <et= 0.00> <mst= 0.00> <sp= 0.00> <cols= 2> <al= 0.00> <DB=-89721.0000> <mlp=10139.0000> <PB=Inf>
<st= 1> <it= 6> <et= 0.00> <mst= 0.00> <sp= 0.00> <cols= 2> <al= 0.00> <DB=-39936.0000> <mlp=10139.0000> <PB=Inf>
<st= 1> <it= 7> <et= 0.00> <mst= 0.00> <sp= 0.00> <cols= 2> <al= 0.00> <DB= -176.0000> <mlp= 124.0000> <PB=Inf>
<st= 1> <it= 8> <et= 0.00> <mst= 0.00> <sp= 0.00> <cols= 2> <al= 0.00> <DB= -6.0000> <mlp= 124.0000> <PB=Inf>
<st= 1> <it= 9> <et= 0.00> <mst= 0.00> <sp= 0.00> <cols= 0> <al= 0.00> <DB= 114.0000> <mlp= 114.0000> <PB=Inf>
***************************************************************************************
**** B&B tree node N°3, parent N°1, depth 1, 2 untreated nodes
**** Local DB = 114.0000, global bounds: [ 114.0000 , Inf ], time = 0.00 sec.
**** Branching constraint: x[1,6]<=0.0
***************************************************************************************
<st= 1> <it= 1> <et= 0.74> <mst= 0.03> <sp= 0.00> <cols= 1> <al= 0.00> <DB= 64.0000> <mlp= 124.0000> <PB=Inf>
<st= 1> <it= 2> <et= 0.74> <mst= 0.00> <sp= 0.00> <cols= 1> <al= 0.00> <DB= 84.0000> <mlp= 124.0000> <PB=Inf>
<st= 1> <it= 3> <et= 0.74> <mst= 0.00> <sp= 0.00> <cols= 0> <al= 0.00> <DB= 114.0000> <mlp= 114.0000> <PB=Inf>
***************************************************************************************
**** B&B tree node N°5, parent N°3, depth 2, 3 untreated nodes
**** Local DB = 114.0000, global bounds: [ 114.0000 , Inf ], time = 0.74 sec.
**** Branching constraint: x[1,5]<=0.0
***************************************************************************************
<st= 1> <it= 1> <et= 0.74> <mst= 0.00> <sp= 0.00> <cols= 1> <al= 0.00> <DB= 64.0000> <mlp= 124.0000> <PB=Inf>
<st= 1> <it= 2> <et= 0.74> <mst= 0.00> <sp= 0.00> <cols= 1> <al= 0.00> <DB= 94.0000> <mlp= 114.0000> <PB=Inf>
<st= 1> <it= 3> <et= 0.74> <mst= 0.00> <sp= 0.00> <cols= 0> <al= 0.00> <DB= 114.0000> <mlp= 114.0000> <PB=Inf>
┌ Warning: No candidate generated. No children will be generated. However, the node is not conquered.
└ @ Coluna.Branching ~/work/Coluna.jl/Coluna.jl/src/Branching/Branching.jl:294
┌ Warning: The solution to the master is not integral and the projection on the original variables is integral.
│ Your reformulation involves subproblems with upper multiplicity greater than 1.
│ Column generation algorithm could not create an integral solution to the master using the column generated.
│ In order to generate columns that can lead to an integral solution, you may have to use a branching scheme that changes the structure of the subproblems.
│ This is not provided by the default implementation of the branching algorithm in the current version of Coluna.
└ @ Coluna.Algorithm ~/work/Coluna.jl/Coluna.jl/src/Algorithm/branching/branchingalgo.jl:104
***************************************************************************************
**** B&B tree node N°4, parent N°3, depth 2, 2 untreated nodes
**** Local DB = 114.0000, global bounds: [ 114.0000 , Inf ], time = 0.74 sec.
**** Branching constraint: x[1,5]>=1.0
***************************************************************************************
<st= 1> <it= 1> <et= 0.74> <mst= 0.00> <sp= 0.00> <cols= 1> <al= 0.00> <DB=-29886.0000> <mlp= 114.0000> <PB=Inf>
<st= 1> <it= 2> <et= 0.74> <mst= 0.00> <sp= 0.00> <cols= 2> <al= 0.00> <DB=-49761.0000> <mlp= 114.0000> <PB=Inf>
<st= 1> <it= 3> <et= 0.74> <mst= 0.00> <sp= 0.00> <cols= 0> <al= 0.00> <DB= 114.0000> <mlp= 114.0000> <PB=Inf>
[ Info: Heuristic Restricted Master IP found improving primal solution with value 114.0
***************************************************************************************
**** B&B tree node N°2, parent N°1, depth 1, 1 untreated node
**** Local DB = 114.0000, global bounds: [ 114.0000 , 114.0000 ], time = 0.74 sec.
**** Branching constraint: x[1,6]>=1.0
***************************************************************************************
──────────────────────────────────────────────────────────────────────────
Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 3.51s / 21.1% 216MiB / 22.7%
Section ncalls time %tot avg alloc %tot avg
──────────────────────────────────────────────────────────────────────────
Coluna 1 743ms 100.0% 743ms 49.1MiB 100.0% 49.1MiB
SolveLpForm 18 33.5ms 4.5% 1.86ms 1.45MiB 3.0% 82.7KiB
──────────────────────────────────────────────────────────────────────────
[ Info: Terminated
[ Info: Primal bound: 114.0
[ Info: Dual bound: 114.0and retrieve the disaggregated solution
for t in T
assignment_patterns = BD.getsolutions(model, t);
for pattern in assignment_patterns
nb_times_pattern_used = BD.value(pattern);
jobs_in_pattern = [];
for j in J
if BD.value(pattern, x[t, j]) ≈ 1
push!(jobs_in_pattern, j);
end
end
println("Pattern of machine type $t used $nb_times_pattern_used times : $jobs_in_pattern");
end
endPattern of machine type 1 used 1.0 times : Any[1, 8]
Pattern of machine type 1 used 1.0 times : Any[2, 4]
Pattern of machine type 1 used 1.0 times : Any[3, 5]
Pattern of machine type 2 used 1.0 times : Any[6, 7]This page was generated using Literate.jl.