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

You 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.01> <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.01> <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.01> <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.01 sec.
**** Branching constraint: x[1,6]<=0.0
***************************************************************************************
  <st= 1> <it=  1> <et= 0.01> <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.01> <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.01> <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.01 sec.
**** Branching constraint: x[1,5]<=0.0
***************************************************************************************
  <st= 1> <it=  1> <et= 0.01> <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.01> <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.01> <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.01 sec.
**** Branching constraint: x[1,5]>=1.0
***************************************************************************************
  <st= 1> <it=  1> <et= 0.01> <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.01> <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.01> <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.01 sec.
**** Branching constraint: x[1,6]>=1.0
***************************************************************************************
 ──────────────────────────────────────────────────────────────────────────
                                  Time                    Allocations
                         ───────────────────────   ────────────────────────
    Tot / % measured:         2.00s /   0.7%            165MiB /  10.3%

 Section         ncalls     time    %tot     avg     alloc    %tot      avg
 ──────────────────────────────────────────────────────────────────────────
 Coluna               1   13.7ms  100.0%  13.7ms   16.9MiB  100.0%  16.9MiB
   SolveLpForm       18   2.40ms   17.5%   133μs    926KiB    5.3%  51.4KiB
 ──────────────────────────────────────────────────────────────────────────
[ Info: Terminated
[ Info: Primal bound: 114.0
[ Info: Dual bound: 114.0

and 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
end

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