Custom Variables and Cuts

Coluna allows users to attach custom data to variables and constraints. This data is useful to store information about the variables or constraints in a custom format much easier to process than extracted information from the formulation (coefficient matrix, bounds, costs, and right-hand side).

In this example, we will show how to attach custom data to variables and constraints and use them to separate non-robust cuts. We will use the Bin Packing problem as an example.

Let us consider a Bin Packing problem with only 3 items such that any pair of items fits into one bin but the 3 items do not. The objective function is to minimize the number of bins being used. Pricing is done by inspection over the 6 combinations of items (3 pairs and 3 singletons). The master LP solution has 1.5 bins at the root node, each 0.5 corresponding to a bin with one of the possible pairs of items.

In this example, we will show you how to use non-robust cuts to improve the master LP solution at the root node. Obviously, Coluna is able to solve this instance by branching on the number of bins but the limit one on the number of nodes prevents it to be solved without cuts.

We define the dependencies:

using JuMP, BlockDecomposition, Coluna, GLPK;

We define the solver.

coluna = JuMP.optimizer_with_attributes(
    Coluna.Optimizer,
    "default_optimizer" => GLPK.Optimizer,
    "params" => Coluna.Params(
        solver = Coluna.Algorithm.TreeSearchAlgorithm(
            conqueralg = Coluna.Algorithm.ColCutGenConquer(
                colgen = Coluna.Algorithm.ColumnGeneration(
                            pricing_prob_solve_alg = Coluna.Algorithm.SolveIpForm(
                                optimizer_id = 1
                            ))
            ),
            maxnumnodes = 1 # we only treat the root node.
        )
    )
);

Let's define the model. Let's $B$ the set of bins and $I$ the set of items. We introduce variable $y_b$ that is equal to 1 if a bin $b$ is used and 0 otherwise. We introduce variable $x_{b,i}$ that is equal to 1 if item $i$ is put in a bin $b$ and 0 otherwise.

model = BlockModel(coluna);

We must assign three items:

I = [1, 2, 3];

And we have three bins:

B = [1, 2, 3];

Each bin is defining a subproblem, we declare our axis:

@axis(axis, collect(B));

We declare subproblem variables y[b]:

@variable(model, y[b in axis], Bin);

And x[b,i]:

@variable(model, x[b in axis, i in I], Bin);

Each item must be assigned to one bin:

@constraint(model, sp[i in I], sum(x[b,i] for b in axis) == 1);

We minimize the number of bins and we declare the decomposition:

@objective(model, Min, sum(y[b] for b in axis))
@dantzig_wolfe_decomposition(model, dec, axis);

Custom data for non-robust cuts

As said previously, at the end of the column generation at the root node, the master LP solution has 1.5 bins. It corresponds to three bins, each of them used 0.5 times containing one pair (1,2), (1, 3), or (2, 3) of items. We are going to introduce the following non-robust cut to make the master LP solution integral:

\[\sum\limits_{s \in S~if~length(s) \geq 2} λ_s \leq 1\]

where :

  • $S$ is the set of possible bin assignments generated by the pricing problem.
  • $length(s)$ the number of items in bin assignment $s \in S$.

This cut means that we cannot have more than one bin with at least two items.

But the problem is that the cut is expressed over the master column and we don't have access to these variables from the JuMP model. To address this problem, Coluna offers a way to compute the coefficient of a column in a constraint by implementing the following method:

Coluna.MathProg.computecoeffFunction
computecoeff(var_custom_data, constr_custom_data) -> Float64

Dispatches on the type of custom data attached to the variable and the constraint to compute the coefficient of the variable in the constraint.

source

We therefore need to attach custom data to the master columns and the non-robust cut to use the method compute_coeff.

For every subproblem solution $s$, we define custom data with the number of items in the bin.

struct MyCustomVarData <: BlockDecomposition.AbstractCustomVarData
    nb_items::Int
end
BlockDecomposition.customvars!(model, MyCustomVarData);

We define custom data for the cut that will contain the minimum number of items in a bin that can be used. The value will be 2 in this example.

struct MyCustomCutData <: BlockDecomposition.AbstractCustomConstrData
    min_items::Int
end
BlockDecomposition.customconstrs!(model, MyCustomCutData);

We implement the computecoeff method for the custom data we defined.

function Coluna.MathProg.computecoeff(
    var_custom_data::MyCustomVarData, constr_custom_data::MyCustomCutData
)
    return (var_custom_data.nb_items >= constr_custom_data.min_items) ? 1.0 : 0.0
end

Pricing callback

We define the pricing callback that will generate the bin with best-reduced cost. Be careful, when using non-robust cuts, you must take into account the contribution of the non-robust cuts to the reduced cost of your solution.

function my_pricing_callback(cbdata)
    # Get the reduced costs of the original variables.
    I = [1, 2, 3]
    b = BlockDecomposition.callback_spid(cbdata, model)

    rc_y = BlockDecomposition.callback_reduced_cost(cbdata, y[b])
    rc_x = [BlockDecomposition.callback_reduced_cost(cbdata, x[b, i]) for i in I]

    # Get the dual values of the custom cuts (to calculate contributions of
    # non-robust cuts to the cost of the solution).
    custduals = Tuple{Int, Float64}[]
    for (_, constr) in Coluna.MathProg.getconstrs(cbdata.form.parent_formulation)
        if typeof(constr.custom_data) == MyCustomCutData
            push!(custduals, (
                constr.custom_data.min_items,
                Coluna.MathProg.getcurincval(cbdata.form.parent_formulation, constr)
            ))
        end
    end

    # Pricing by inspection.
    sols = [[1], [2], [3], [1, 2], [1, 3], [2, 3]]
    best_s = Int[]
    best_rc = Inf
    for s in sols
        rc_s = rc_y + sum(rc_x[i] for i in s) # reduced cost of the subproblem variables
        if !isempty(custduals)
            # contribution of the non-robust cuts
            rc_s -= sum((length(s) >= minits) ? dual : 0.0 for (minits, dual) in custduals)
        end
        if rc_s < best_rc
            best_rc = rc_s
            best_s = s
        end
    end
    @show best_s
    # build the best one and submit
    solcost = best_rc
    solvars = JuMP.VariableRef[]
    solvarvals = Float64[]
    for i in best_s
        push!(solvars, x[b, i])
        push!(solvarvals, 1.0)
    end
    push!(solvars, y[b])
    push!(solvarvals, 1.0)
    # submit the solution
    MOI.submit(
        model, BlockDecomposition.PricingSolution(cbdata),
        solcost,
        solvars,
        solvarvals,
        MyCustomVarData(length(best_s)) # attach a custom data to the column
    )

    MOI.submit(model, BlockDecomposition.PricingDualBound(cbdata), solcost)
    return
end
my_pricing_callback (generic function with 1 method)

The pricing callback is done, we define it as the solver of our pricing problem.

subproblems = BlockDecomposition.getsubproblems(dec)
BlockDecomposition.specify!.(
    subproblems,
    solver = my_pricing_callback
);

Non-robust cut separation callback.

We now define the cut separation callback for our non-robust cut. This is the same callback as the one used for robust cuts. There is just one slight difference when you submit the non-robust cut. Since cuts are expressed over the master variables and these variables are inaccessible from the JuMP model, you'll submit a constraint with an empty left-hand side and you'll leave Coluna populate the left-hand side with the values returned by Coluna.MathProg.computecoeff.

So let's define the callback. Basically, if the solution uses more than one bin with two items, The cut is added to the model.

function custom_cut_sep(cbdata)
    # Compute the constraint violation by iterating over the master solution.
    viol = -1.0
    for (varid, varval) in cbdata.orig_sol
        var = Coluna.MathProg.getvar(cbdata.form, varid)
        if !isnothing(var.custom_data)
            if var.custom_data.nb_items >= 2
                viol += varval
            end
        end
    end
    # Add the cut (at most one variable with 2 or more of the 3 items) if violated.
    if viol > 0.001
        MOI.submit(
            model, MOI.UserCut(cbdata),
            JuMP.ScalarConstraint(
                JuMP.AffExpr(0.0), # We cannot express the left-hand side so we push 0.
                MOI.LessThan(1.0)
            ),
            MyCustomCutData(2) # Cut custom data.
        )
    end
    return
end

MOI.set(model, MOI.UserCutCallback(), custom_cut_sep)
JuMP.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.
***************************************************************************************
best_s = [1, 2]
best_s = [1, 2]
best_s = [1, 2]
  <st= 1> <it=  1> <et= 2.09> <mst= 0.01> <sp= 0.00> <cols= 3> <al= 0.00> <DB=-29997.0000> <mlp=60000.0000> <PB=Inf>
best_s = [1, 3]
best_s = [1, 3]
best_s = [1, 3]
  <st= 1> <it=  2> <et= 2.12> <mst= 0.03> <sp= 0.00> <cols= 3> <al= 0.00> <DB=-29999.0000> <mlp=30001.0000> <PB=Inf>
best_s = [2, 3]
best_s = [2, 3]
best_s = [2, 3]
  <st= 1> <it=  3> <et= 2.12> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB=-10001.0000> <mlp=20002.0000> <PB=Inf>
best_s = [1]
best_s = [1]
best_s = [1]
  <st= 1> <it=  4> <et= 2.12> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB=    3.0000> <mlp=15001.5000> <PB=Inf>
best_s = [3]
best_s = [3]
best_s = [3]
  <st= 1> <it=  5> <et= 2.12> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB=-19995.0000> <mlp=10002.0000> <PB=Inf>
best_s = [2]
best_s = [2]
best_s = [2]
  <st= 1> <it=  6> <et= 2.12> <mst= 0.00> <sp= 0.00> <cols= 3> <al= 0.00> <DB=-19995.0000> <mlp=10002.0000> <PB=Inf>
Cut separation callback adds 0 new essential cuts and 0 new facultative cuts.
[ Info: Improving primal solution with value 3.0 is found during column generation
best_s = [1]
best_s = [1]
best_s = [1]
  <st= 1> <it=  7> <et= 2.12> <mst= 0.00> <sp= 0.00> <cols= 0> <al= 0.00> <DB=    3.0000> <mlp=    3.0000> <PB=3.0000>
───────────────────────────────────────────────────────────────────────────────
                                      Time                    Allocations
                             ───────────────────────   ────────────────────────
      Tot / % measured:           5.68s /  37.3%           2.21GiB /   5.1%

Section              ncalls     time    %tot     avg     alloc    %tot      avg
───────────────────────────────────────────────────────────────────────────────
Coluna                    1    2.12s  100.0%   2.12s    117MiB  100.0%   117MiB
  Pricing callback       21    352ms   16.6%  16.8ms   31.0MiB   26.6%  1.48MiB
  SolveLpForm             7   38.3ms    1.8%  5.48ms   1.85MiB    1.6%   271KiB
───────────────────────────────────────────────────────────────────────────────
[ Info: Terminated
[ Info: Primal bound: 3.0
[ Info: Dual bound: 3.0

We see on the output that the algorithm has converged a first time before a cut is added. Coluna then starts a new iteration taking into account the cut. We notice here an improvement of the value of the dual bound: before the cut, we converge towards 1.5. After the cut, we reach 2.0.


This page was generated using Literate.jl.