Modeling Packing Problems with Set Decisions¶

The Bin Packing Problem (BPP)¶

In the Bin Packing Problem (BPP), a number of items with known weights must be assigned to bins with uniform capacity. Each item must be placed inside exactly one bin, and the total weight of items inside each bin must not exceed its capacity. The objective is to minimize the number of bins used.

While a classic MILP model using Boolean decisions exists for the Bin Packing Problem, we show how the set decisions in Hexaly’ modeling formalism enable users to write a more straightforward and more compact model. The difference in the complexity of the models makes the set-based model all the more useful that it enables the solver to scale to much larger instances, and deliver much better performance than classic MILP solvers.

We start by defining the decisions. The classic MILP model uses Boolean decisions (x[i][k] is true if item i is assigned to bin k):

Instead of Booleans, the Hexaly model for the Bin Packing Problem (BPP) uses set decision variables. Mathematically, a set decision variable of domain size n represents a subset of {0, ..., n-1}. It is a set in the usual mathematical sense: each element is unique in the set, and the elements inside a set are unordered. In the Bin Packing Problem, each set represents a bin: the elements inside bins[k] are the indices of the items assigned to bin k:

bins[k in 0...nbMaxBins] <- set(nbItems);

bins = [model.set(nb_items) for _ in range(nb_max_bins)]

for (int k = 0; k < nbMaxBins; ++k) {
    bins[k] = model.setVar(nbItems);
}

for (int k = 0; k < nbMaxBins; ++k)
    bins[k] = model.Set(nbItems);

for (int k = 0; k < nbMaxBins; ++k) {
    bins[k] = model.setVar(nbItems);
}

The model then contains O(n) decisions, instead of O(n²) for the Boolean model.

After defining the decision variables, we must ensure that each item is assigned to exactly one bin. The Boolean model needs n constraints to do so:

In the set-based model, only one constraint is needed thanks to the partition operator:

constraint partition[k in 0...nbMaxBins](bins[k]);

model.constraint(model.partition(bins))

model.constraint(model.partition(bins.begin(), bins.end()));

model.Constraint(model.Partition(bins));

model.constraint(model.partition(bins));

The partition operator, when applied to a series of sets of domain size n, returns 1 if every element between 0 and n-1 is present in exactly one of the sets. Therefore, this constraint ensures that the index of each item is present in exactly one of the sets, meaning that each item is assigned to exactly one of the bins.

We then need to enforce the capacity constraints: the total weight of items in each bin must not exceed the bin’s capacity. In the Boolean model, we need to compute n sums with n terms each:

In the set-based model, we can reduce the complexity of the weight computation by using the variadic form of the sum operator along with a lambda function returning the weight of an item. In the following code, binWeights is not a decision variable but an intermediate expression: its value is directly computed from the values of the decisions and the input data:

binWeights[k] <- sum(bins[k], i => itemWeights[i]);
constraint binWeights[k] <= binCapacity;

weight_lambda = model.lambda_function(lambda i: weights[i])
bin_weights = [model.sum(b, weight_lambda) for b in bins]
for w in bin_weights:
    model.constraint(w <= bin_capacity)

HxExpression weightLambda = model.createLambdaFunction([&](HxExpression i) { return weights[i]; });
for (int k = 0; k < nbMaxBins; ++k) {
    binWeights[k] = model.sum(bins[k], weightLambda);
    model.constraint(binWeights[k] <= binCapacity);
}

HxExpression weightLambda = model.LambdaFunction(i => weights[i]);
for (int k = 0; k < nbMaxBins; ++k)
{
    binWeights[k] = model.Sum(bins[k], weightLambda);
    model.Constraint(binWeights[k] <= binCapacity);
}

HxExpression weightLambda = model.lambdaFunction(i -> model.at(weights, i));
for (int k = 0; k < nbMaxBins; ++k) {
    binWeights[k] = model.sum(bins[k], weightLambda);
    model.constraint(model.leq(binWeights[k], binCapacity));
}

This expression loops over all the elements in bins[k] to sum the associated weights. Note that the number of terms in this sum varies during the search (it is equal to the number of items assigned to bin k). This feature of Hexaly Optimizer’s modeling formalism is particularly powerful, since being able to apply n-ary operators to a variable number of terms can make the models much more compact, and lead to better performance overall on large instances.

Indeed, since each element belongs to exactly one bin, each element is only involved in the computation of one sum. The total complexity for computing the weights of all the bins is then O(n) for the set-based model, as opposed to O(n²) for the Boolean model.

To compute the objective function, we need to check whether each bin is actually used or not. For that, the Boolean model uses n sums again:

In the set-based model, we simply use the count operator, which returns the number of elements in a set:

binsUsed[k] <- (count(bins[k]) > 0);

bins_used = [model.count(b) > 0 for b in bins]

binsUsed[k] = model.count(bins[k]) > 0;

binsUsed[k] = model.Count(bins[k]) > 0;

binsUsed[k] = model.gt(model.count(bins[k]), 0);

The complete Hexaly model for the Bin Packing Problem (BPP) is then:

function model() {
    bins[k in 0...nbMaxBins] <- set(nbItems);
    constraint partition[k in 0...nbMaxBins](bins[k]);

    for [k in 0...nbMaxBins] {
        binWeights[k] <- sum(bins[k], i => itemWeights[i]);
        constraint binWeights[k] <= binCapacity;

        binsUsed[k] <- (count(bins[k]) > 0);
    }
    totalBinsUsed <- sum[k in 0...nbMaxBins](binsUsed[k]);
    minimize totalBinsUsed;
}

bins = [model.set(nb_items) for _ in range(nb_max_bins)]
model.constraint(model.partition(bins))

weights = model.array(weights_data)
weight_lambda = model.lambda_function(lambda i: weights[i])

bin_weights = [model.sum(b, weight_lambda) for b in bins]
for w in bin_weights:
    model.constraint(w <= bin_capacity)

bins_used = [model.count(b) > 0 for b in bins]

total_bins_used = model.sum(bins_used)
model.minimize(total_bins_used)

for (int k = 0; k < nbMaxBins; ++k) {
    bins[k] = model.setVar(nbItems);
}
model.constraint(model.partition(bins.begin(), bins.end()));

HxExpression weights = model.array(weightsData.begin(), weightsData.end());
HxExpression weightLambda = model.createLambdaFunction([&](HxExpression i) { return weights[i]; });

for (int k = 0; k < nbMaxBins; ++k) {
    binWeights[k] = model.sum(bins[k], weightLambda);
    model.constraint(binWeights[k] <= binCapacity);

    binsUsed[k] = model.count(bins[k]) > 0;
}

totalBinsUsed = model.sum(binsUsed.begin(), binsUsed.end());
model.minimize(totalBinsUsed);

for (int k = 0; k < nbMaxBins; ++k)
    bins[k] = model.Set(nbItems);
model.Constraint(model.Partition(bins));

HxExpression weights = model.Array(weightsData);
HxExpression weightLambda = model.LambdaFunction(i => weights[i]);

for (int k = 0; k < nbMaxBins; ++k)
{
    binWeights[k] = model.Sum(bins[k], weightLambda);
    model.Constraint(binWeights[k] <= binCapacity);

    binsUsed[k] = model.Count(bins[k]) > 0;
}

totalBinsUsed = model.Sum(binsUsed);
model.Minimize(totalBinsUsed);

for (int k = 0; k < nbMaxBins; ++k) {
    bins[k] = model.setVar(nbItems);
}
model.constraint(model.partition(bins));

HxExpression weights = model.array(weightsData);
HxExpression weightLambda = model.lambdaFunction(i -> model.at(weights, i));

for (int k = 0; k < nbMaxBins; ++k) {
    binWeights[k] = model.sum(bins[k], weightLambda);
    model.constraint(model.leq(binWeights[k], binCapacity));

    binsUsed[k] = model.gt(model.count(bins[k]), 0);
}

totalBinsUsed = model.sum(binsUsed);
model.minimize(totalBinsUsed);

The whole code for solving the Bin Packing Problem (BPP) with Hexaly, along with a few instances, is available as a code template, with an associated BPP benchmark.

The Bin Packing Problem with Conflicts (BPPC)¶

The Bin Packing Problem with Conflicts (BPPC) is an extension of the Bin Packing Problem (BPP), in which certain pairs of items conflict and, therefore, cannot be placed in the same bin.

To write the conflict constraints, we first need to retrieve the index of the bin containing each item, using the find operator:

for [i in 0...nbItems] {
    binForItem[i] <- find(bins, i);
}

bin_for_item = [model.find(bins_array, i) for i in range(nb_items)]

for (int i = 0; i < nbItems; ++i) {
    binForItem[i] = model.find(binsArray, i);
}

for (int i = 0; i < nbItems; ++i) {
    binForItem[i] = model.Find(binsArray, i);
}

for (int i = 0; i < nbItems; ++i) {
    binForItem[i] = model.find(binsArray, i);
}

In the expression above, the find operator returns the position of the first set containing element i.

To enforce the conflict constraints, we could simply impose that for any pair of items (i, j) that conflict, binForItem[i] and binForItem[j] must be different. However, such a formulation would result in O(n²) constraints, which is not ideal for large-scale instances.

Instead, using the intersection and count operators, we ensure that the intersection between an item’s chosen bin and its set of conflicting items is empty:

for [i in 0...nbItems] {
    itemsIntersection  <- intersection(bins[binForItem[i]], forbiddenItems[i]);
    constraint count(itemsIntersection) == 0;
}

for i in range(nb_items):
    items_intersection = model.intersection(forbidden_items_array[i], bins_array[bin_for_item[i]])
    model.constraint(model.count(items_intersection) == 0)

for (int i = 0; i < nbItems; ++i) {
    HxExpression itemsIntersection = model.intersection(binsArray[binForItem[i]], forbiddenItemsArray[i]);
    model.constraint(model.count(itemsIntersection) == 0);
}

for (int i = 0; i < nbItems; ++i)
{
    HxExpression itemsIntersection = model.Intersection(binsArray[binForItem[i]], forbiddenItemsArray[i]);
    model.Constraint(model.Count(itemsIntersection) == 0);
}

for (int i = 0; i < nbItems; ++i) {
    HxExpression itemsIntersection = model.intersection(model.at(binsArray, binForItem[i]),
            model.at(forbiddenItemsArray, i));
    model.constraint(model.eq(model.count(itemsIntersection), 0));
}

The formulation above is much more compact, with only O(n) constraints, and helps the solver scale to much larger instances.

The complete Hexaly model for the Bin Packing Problem with Conflicts (BPPC), along with a few instances, is available as a code template, with an associated BPPC benchmark.

The Constrained Pit Limit Problem (CPIT)¶

The Constrained Pit Limit Problem (CPIT) is an extension of the Bin Packing Problem (BPP) with precedence constraints. In the Constrained Pit Limit Problem (CPIT), we consider a set of blocks that can be extracted from a mining pit during several time periods. Extracting a block requires a certain amount of resources. For each type of resource, the sum of resources consumed over a period cannot exceed a certain limit. Each block must be extracted during the same period or during a later period than all its predecessors. The profit of a block’s extraction depends on the block and the period of extraction. Finally, the objective function is to maximize profits made from block extractions.

Similarly to the Bin Packing Problem (BPP), the Constrained Pit Limit Problem (CPIT) uses set decision variables, representing the blocks extracted during each period (including a dummy period at the end, containing all unextracted blocks). To ensure that all blocks belong to exactly one period, the sets are constrained to form a partition:

periods[t in 0...nbPeriods+1] <- set(nbBlocks);
constraint partition(periods);

period_vars = [model.set(nb_blocks) for t in range(nb_periods+1)]
periods = model.array(period_vars)
model.constraint(model.partition(periods))

HxExpression periods = model.array();
for (int t=0; t<nbPeriods+1; t++){
    periodVars[t] = model.setVar(nbBlocks);
    periods.addOperand(periodVars[t]);
}
model.constraint(model.partition(periods));

HxExpression periods = model.Array();
for (int t = 0; t < nbPeriods + 1; t++)
{
    periodVars[t] = model.Set(nbBlocks);
    periods.AddOperand(periodVars[t]);
}
model.Constraint(model.Partition(periods));

HxExpression periods = model.array();
for (int t = 0; t < nbPeriods + 1; t++){
    periodVars[t] = model.setVar(nbBlocks);
    periods.addOperand(periodVars[t]);
}
model.constraint(model.partition(periods));

The resource constraints are very similar to the capacity constraints in the Bin Packing Problem (BPP). We compute the resource consumption at each period for each type of resource using the variadic form of the sum operator along with a lambda function returning the resource use of a block:

sumResourceConsumption[r in 0...nbResources][t in 0...nbPeriods]
        <- sum(periods[t], b => resourceUse[r][b]);

consumption_lambda = [model.lambda_function(lambda b: resource_use[r][b]) \
        for r in range(nb_resources)]
consumption_per_period = [[model.sum(period_vars[t],consumption_lambda[r]) \
        for t in range(nb_periods)]  for r in range(nb_resources)]

for (int r = 0; r < nbResources; r++){
    consumptionLambda.push_back(model.lambdaFunction([&](HxExpression b)
            { return resourceUse[r][b]; }));
    for (int t=0; t<nbPeriods; t++){
        sumResourceConsumption[r][t] = model.sum(periodVars[t], consumptionLambda[r]);
    }
}

for (int r = 0; r < nbResources; r++)
{
    consumptionLambda[r] = model.LambdaFunction(b =>
    {
        return resourceUse[r][b];
    });
    for (int t = 0; t < nbPeriods; t++)
    {
        sumResourceConsumption[r][t] = model.Sum(periodVars[t], consumptionLambda[r]);
    }
}

for (int r = 0; r < nbResources; r++){
    HxExpression resourceUseR = resourceUse[r];
    consumptionLambda[r] = model.lambdaFunction(b -> model.at(resourceUseR,b));
    for (int t = 0; t < nbPeriods; t++){
        sumResourceConsumption[r][t] = model.sum(periodVars[t], consumptionLambda[r]);
    }
}

This expression loops over all the blocks in periods[t] to sum the associated resource use. Note that the number of terms in this sum varies during the search (it is equal to the number of blocks extracted at period t). We can then constrain it to be above the minimum resource consumption and below the maximum resource consumption.

To write the precedence constraints between blocks, we need to retrieve each block’s extraction period. Like in the Bin Packing Problem with Conflicts (BPPC), we use the find operator:

blockExtractionPeriod[b in 0...nbBlocks] <- find(periods, b);

block_extraction_period = [model.find(periods, b) for b in range(nb_blocks)]

for (int b=0; b<nbBlocks; b++){
    blockExtractionPeriod[b] = model.find(periods,b);
}

for (int b = 0; b < nbBlocks; b++)
    blockExtractionPeriod[b] = model.Find(periods, b);

for (int b = 0; b < nbBlocks; b++){
    blockExtractionPeriod[b] = model.find(periods, b);
}

For each block and each predecessor of this block, we ensure that the block’s extraction period has an index greater than or equal to that of the predecessor’s extraction period:

for [b in 0...nbBlocks]{
    for [pred in predecessors[b]]{
        constraint blockExtractionPeriod[b] >= blockExtractionPeriod[pred];
    }
}

for b in range(nb_blocks):
    for pred in precedence[b]:
        model.constraint(block_extraction_period[b] >= block_extraction_period[pred])

for (int b = 0; b < nbBlocks; b++){
    for (int pred: precedence[b]){
        model.constraint(blockExtractionPeriod[b] >= blockExtractionPeriod[pred]);
    }
}

for (int b = 0; b < nbBlocks; b++)
{
    foreach (int pred in precedence[b])
    {
        model.Constraint(blockExtractionPeriod[b] >= blockExtractionPeriod[pred]);
    }
}

for (int b = 0; b < nbBlocks; b++){
    ArrayList<Integer> predecessors_b = precedence.get(b);
    for (int pred : predecessors_b){
        model.constraint(model.leq(blockExtractionPeriod[pred],blockExtractionPeriod[b]));
    }
}

An alternative, more compact formulation replaces these constraints with one constraint per block with predecessors: the block’s extraction period must be greater than or equal to the maximum extraction period of its predecessors. This enforces the same precedence relations with fewer constraints and can improve performance:

for [b in 0...nbBlocks : predecessors[b].count() > 0] {
    constraint blockExtractionPeriod[b]
            >= max[pred in predecessors[b]](blockExtractionPeriod[pred]);
}

for b in range(nb_blocks):
    if precedence.get(b):
        latest_predecessor_period = model.max(
            block_extraction_period[pred] for pred in precedence[b])
        model.constraint(block_extraction_period[b] >= latest_predecessor_period)

for (int b = 0; b < nbBlocks; ++b) {
    if (!precedence[b].empty()) {
        HxExpression latestPredecessorPeriod = model.max();
        for (int pred : precedence[b]) {
            latestPredecessorPeriod.addOperand(blockExtractionPeriod[pred]);
        }
        model.constraint(blockExtractionPeriod[b] >= latestPredecessorPeriod);
    }
}

for (int b = 0; b < nbBlocks; ++b)
{
    if (precedence[b] != null && precedence[b].Count > 0)
    {
        HxExpression latestPredecessorPeriod = model.Max();
        foreach (int pred in precedence[b])
        {
            latestPredecessorPeriod.AddOperand(blockExtractionPeriod[pred]);
        }
        model.Constraint(blockExtractionPeriod[b] >= latestPredecessorPeriod);
    }
}

for (int b = 0; b < nbBlocks; ++b) {
    ArrayList<Integer> predecessors = precedence.get(b);
    if (!predecessors.isEmpty()) {
        HxExpression latestPredecessorPeriod = model.max();
        for (int pred : predecessors) {
            latestPredecessorPeriod.addOperand(blockExtractionPeriod[pred]);
        }
        model.constraint(model.geq(blockExtractionPeriod[b], latestPredecessorPeriod));
    }
}

The complete Hexaly model for the Constrained Pit Limit Problem (CPIT), along with a few instances, is available as a code template.

Further reading: discover our code templates on packing problems¶