Column generation with the Generalized Assignment Problem
This quick start guide introduces the main features of Coluna through the example of the Generalized Assignment Problem.
Classic model solved with MIP solver
Consider a set of machines M and a set of jobs J. A machine $m$ has a resource capacity $Q_m$ . A job $j$ assigned to a machine $m$ has a cost $c_{mj}$ and consumes $w_{mj}$ resource units of the machine $m$. The goal is to minimize the sum of job costs while assigning each job to a machine and not exceeding the capacity of each machine.
Let $x_{mj}$ equal to one if job $j$ is assigned to machine $m$; $0$ otherwise. The problem has the original formulation:
\[\begin{alignedat}{4} \text{[GAP]} \equiv \min \mathrlap{\sum_{m \in M}\sum_{j \in J} c_{mj} x_{mj}} \\ \text{s.t.} && \sum_{m \in M} x_{mj} &= 1 \quad& j \in J \\ && \sum_{j \in J} w_{mj} x_{mj} &\leq Q_m \quad \quad& m \in M \\ && x_{mj} &\in \{0,1\} &m \in M,\; j \in J \end{alignedat}\]
Let us consider the following instance.
M = 1:3;
J = 1:15;
c = [12.7 22.5 8.9 20.8 13.6 12.4 24.8 19.1 11.5 17.4 24.7 6.8 21.7 14.3 10.5; 19.1 24.8 24.4 23.6 16.1 20.6 15.0 9.5 7.9 11.3 22.6 8.0 21.5 14.7 23.2; 18.6 14.1 22.7 9.9 24.2 24.5 20.8 12.9 17.7 11.9 18.7 10.1 9.1 8.9 7.7; 13.1 16.2 16.8 16.7 9.0 16.9 17.9 12.1 17.5 22.0 19.9 14.6 18.2 19.6 24.2];
w = [61 70 57 82 51 74 98 64 86 80 69 79 60 76 78; 50 57 61 83 81 79 63 99 82 59 83 91 59 99 91;91 81 66 63 59 81 87 90 65 55 57 68 92 91 86; 62 79 73 60 75 66 68 99 69 60 56 100 67 68 54];
Q = [1020 1460 1530];We write the model with JuMP, a domain-specific modeling language for mathematical optimization embedded in Julia. We optimize with GLPK. If you are not familiar with the JuMP package, you may want to check its documentation.
using JuMP, GLPK;A JuMP model for the original formulation is:
model = Model(GLPK.Optimizer)
@variable(model, x[m in M, j in J], Bin);
@constraint(model, cov[j in J], sum(x[m, j] for m in M) >= 1);
@constraint(model, knp[m in M], sum(w[m, j] * x[m, j] for j in J) <= Q[m]);
@objective(model, Min, sum(c[m, j] * x[m, j] for m in M, j in J));We optimize the instance and retrieve the objective value.
optimize!(model);
objective_value(model)166.5Try column generation easily with Coluna and BlockDecomposition
This model has a block structure: each knapsack constraint defines an independent block and the set-partitioning constraints couple these independent blocks. By applying the Dantzig-Wolfe reformulation, each knapsack constraint forms a tractable subproblem and the set-partitioning constraints are handled in a master problem.
To write the model, you need JuMP and BlockDecomposition. The latter is an extension built on top of JuMP to model Dantzig-Wolfe and Benders decompositions. You will find more documentation about BlockDecomposition in the Decomposition & reformulation To optimize the problem, you need Coluna and a Julia package that provides a MIP solver such as GLPK.
Since we have already loaded JuMP and GLPK, we just need:
using BlockDecomposition, Coluna;Next, you instantiate the solver and define the algorithm that you use to optimize the problem. In this case, the algorithm is a classic branch-and-price provided by Coluna.
coluna = optimizer_with_attributes(
Coluna.Optimizer,
"params" => Coluna.Params(
solver = Coluna.Algorithm.TreeSearchAlgorithm() # default branch-cut-and-price
),
"default_optimizer" => GLPK.Optimizer # GLPK for the master & the subproblems
);In BlockDecomposition, an axis is an index set of subproblems. Let M_axis be the index set of machines; it defines an axis along which we can implement the desired decomposition.
@axis(M_axis, M);In this example, the axis M_axis defines one knapsack subproblem for each machine. For instance, the first machine index is 1 and is of type BlockDecomposition.AxisId:
M_axis[1]
typeof(M_axis[1])BlockDecomposition.AxisId{:M_axis, Int64}Jobs are not involved in the decomposition, set J of jobs thus stays as a classic range.
The model takes the form:
model = BlockModel(coluna);You can write BlockModel(coluna; direct_model = true) to pass names of variables and constraints to Coluna.
@variable(model, x[m in M_axis, j in J], Bin);
@constraint(model, cov[j in J], sum(x[m, j] for m in M_axis) >= 1);
@constraint(model, knp[m in M_axis], sum(w[m, j] * x[m, j] for j in J) <= Q[m]);
@objective(model, Min, sum(c[m, j] * x[m, j] for m in M_axis, j in J));This is the same model as above except that we use a BlockModel instead of a Model and M_axis as the set of machines instead of M. Therefore, BlockDecomposition will know which variables and constraints are involved in subproblems because one of their indices is a BlockDecomposition.AxisId.
You then apply a Dantzig-Wolfe decomposition along M_axis:
@dantzig_wolfe_decomposition(model, decomposition, M_axis)where decomposition is a variable that contains information about the decomposition.
decompositionRoot - Annotation(BlockDecomposition.Master, BlockDecomposition.DantzigWolfe, lm = 1.0, um = 1.0, bp = 1.0, id = 2) with 3 subproblems :
2 => Annotation(BlockDecomposition.DwPricingSp, BlockDecomposition.DantzigWolfe, lm = 1.0, um = 1.0, bp = 1.0, id = 4)
3 => Annotation(BlockDecomposition.DwPricingSp, BlockDecomposition.DantzigWolfe, lm = 1.0, um = 1.0, bp = 1.0, id = 5)
1 => Annotation(BlockDecomposition.DwPricingSp, BlockDecomposition.DantzigWolfe, lm = 1.0, um = 1.0, bp = 1.0, id = 3)
Once the decomposition is defined, you can retrieve the master and the subproblems to give additional information to the solver.
master = getmaster(decomposition)
subproblems = getsubproblems(decomposition)3-element Vector{BlockDecomposition.SubproblemForm}:
Subproblem formulation for M_axis = 1 contains : 1.0 <= multiplicity <= 1.0
Subproblem formulation for M_axis = 2 contains : 1.0 <= multiplicity <= 1.0
Subproblem formulation for M_axis = 3 contains : 1.0 <= multiplicity <= 1.0
The multiplicity of a subproblem is the number of times that the same independent block shaped by the subproblem appears in the model. This multiplicity also specifies the number of solutions to the subproblem that can appear in the solution to the original problem.
In this GAP instance, the upper multiplicity is $1$ because every subproblem is different, i.e., every machine is different and used at most once.
The lower multiplicity is $0$ because a machine may stay unused. The multiplicity specifications take the form:
specify!.(subproblems, lower_multiplicity = 0, upper_multiplicity = 1)
getsubproblems(decomposition)3-element Vector{BlockDecomposition.SubproblemForm}:
Subproblem formulation for M_axis = 1 contains : 0.0 <= multiplicity <= 1.0
Subproblem formulation for M_axis = 2 contains : 0.0 <= multiplicity <= 1.0
Subproblem formulation for M_axis = 3 contains : 0.0 <= multiplicity <= 1.0
The model is now fully defined. To solve it, you need to call:
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= 3> <al= 0.00> <DB=-199514.2000> <mlp=100000.0000> <PB=Inf>
<st= 1> <it= 2> <et= 0.00> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB=-99999.8000> <mlp=50148.8000> <PB=Inf>
[ Info: Improving primal solution with value 220.2 is found during column generation
<st= 1> <it= 3> <et= 0.00> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= -453.9000> <mlp= 220.2000> <PB=220.2000>
<st= 1> <it= 4> <et= 0.00> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= -411.4000> <mlp= 220.2000> <PB=220.2000>
<st= 1> <it= 5> <et= 0.00> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= -350.0000> <mlp= 220.2000> <PB=220.2000>
<st= 1> <it= 6> <et= 0.00> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= -314.8000> <mlp= 220.2000> <PB=220.2000>
<st= 1> <it= 7> <et= 0.00> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= -299.3000> <mlp= 220.2000> <PB=220.2000>
<st= 1> <it= 8> <et= 0.00> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= -285.0000> <mlp= 220.2000> <PB=220.2000>
<st= 1> <it= 9> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= -155.5667> <mlp= 220.2000> <PB=220.2000>
<st= 1> <it= 10> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= -139.9400> <mlp= 220.2000> <PB=220.2000>
<st= 1> <it= 11> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= -63.5455> <mlp= 220.2000> <PB=220.2000>
<st= 1> <it= 12> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= -51.8000> <mlp= 220.2000> <PB=220.2000>
<st= 1> <it= 13> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 61.1900> <mlp= 208.7500> <PB=220.2000>
<st= 1> <it= 14> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= -13.0964> <mlp= 207.9821> <PB=220.2000>
<st= 1> <it= 15> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 57.6200> <mlp= 206.7200> <PB=220.2000>
<st= 1> <it= 16> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 149.3588> <mlp= 204.8441> <PB=220.2000>
<st= 1> <it= 17> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 95.3995> <mlp= 202.4207> <PB=220.2000>
<st= 1> <it= 18> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 107.0038> <mlp= 196.0509> <PB=220.2000>
<st= 1> <it= 19> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 78.5000> <mlp= 188.6667> <PB=220.2000>
<st= 1> <it= 20> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 110.5091> <mlp= 186.5545> <PB=220.2000>
<st= 1> <it= 21> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 150.2941> <mlp= 185.3540> <PB=220.2000>
[ Info: Improving primal solution with value 180.6 is found during column generation
<st= 1> <it= 22> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 113.6444> <mlp= 180.6000> <PB=180.6000>
<st= 1> <it= 23> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 121.1800> <mlp= 180.6000> <PB=180.6000>
<st= 1> <it= 24> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 121.5385> <mlp= 177.6385> <PB=180.6000>
<st= 1> <it= 25> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 135.5560> <mlp= 176.6960> <PB=180.6000>
<st= 1> <it= 26> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 124.8000> <mlp= 175.4000> <PB=180.6000>
<st= 1> <it= 27> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 140.3645> <mlp= 175.3548> <PB=180.6000>
<st= 1> <it= 28> <et= 0.01> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 140.6211> <mlp= 174.9184> <PB=180.6000>
<st= 1> <it= 29> <et= 0.02> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 131.2406> <mlp= 173.4125> <PB=180.6000>
<st= 1> <it= 30> <et= 0.02> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 146.2667> <mlp= 173.0521> <PB=180.6000>
<st= 1> <it= 31> <et= 0.02> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 152.7125> <mlp= 173.0521> <PB=180.6000>
<st= 1> <it= 32> <et= 0.02> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 155.1333> <mlp= 172.6444> <PB=180.6000>
<st= 1> <it= 33> <et= 0.02> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 156.4722> <mlp= 172.6444> <PB=180.6000>
[ Info: Improving primal solution with value 171.10000000000002 is found during column generation
<st= 1> <it= 34> <et= 0.02> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 159.0324> <mlp= 171.1000> <PB=171.1000>
[ Info: Improving primal solution with value 167.70000000000002 is found during column generation
<st= 1> <it= 35> <et= 0.02> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB= 163.5000> <mlp= 167.7000> <PB=167.7000>
<st= 1> <it= 36> <et= 0.02> <mst= 0.00> <sp= 0.00> <cols= 1> <al= 0.00> <DB= 164.4000> <mlp= 167.7000> <PB=167.7000>
<st= 1> <it= 37> <et= 0.02> <mst= 0.00> <sp= 0.00> <cols= 1> <al= 0.00> <DB= 165.9000> <mlp= 167.7000> <PB=167.7000>
<st= 1> <it= 38> <et= 0.02> <mst= 0.00> <sp= 0.00> <cols= 1> <al= 0.00> <DB= 166.5000> <mlp= 167.7000> <PB=167.7000>
<st= 1> <it= 39> <et= 0.02> <mst= 0.00> <sp= 0.00> <cols= 2> <al= 0.00> <DB= 165.1000> <mlp= 167.7000> <PB=167.7000>
<st= 1> <it= 40> <et= 0.02> <mst= 0.00> <sp= 0.00> <cols= 1> <al= 0.00> <DB= 164.5000> <mlp= 167.7000> <PB=167.7000>
<st= 1> <it= 41> <et= 0.02> <mst= 0.00> <sp= 0.00> <cols= 1> <al= 0.00> <DB= 166.3286> <mlp= 167.7000> <PB=167.7000>
<st= 1> <it= 42> <et= 0.02> <mst= 0.00> <sp= 0.00> <cols= 1> <al= 0.00> <DB= 166.5000> <mlp= 167.7000> <PB=167.7000>
[ Info: Improving primal solution with value 166.5 is found during column generation
<st= 1> <it= 43> <et= 0.02> <mst= 0.00> <sp= 0.00> <cols= 0> <al= 0.00> <DB= 166.5000> <mlp= 166.5000> <PB=166.5000>
───────────────────────────────────────────────────────────────────────────────
Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 1.86s / 28.8% 240MiB / 46.8%
Section ncalls time %tot avg alloc %tot avg
───────────────────────────────────────────────────────────────────────────────
Coluna 2 536ms 100.0% 268ms 112MiB 100.0% 56.1MiB
Pricing callback 143 487ms 91.0% 3.41ms 69.7MiB 62.1% 499KiB
SolveLpForm 91 21.6ms 4.0% 237μs 10.7MiB 9.5% 120KiB
───────────────────────────────────────────────────────────────────────────────
[ Info: Terminated
[ Info: Primal bound: 166.5
[ Info: Dual bound: 166.5You can find more information about the output of the column generation algorithm ColumnGeneration.
Finally, you can retrieve the solution to the original formulation with JuMP methods. For example, if we want to know if job 3 is assigned to machine 1:
value(x[1,3])1.0This page was generated using Literate.jl.