Stochastic Packing¶

Principles learned¶

  • Use list decision variables

  • Use a lambda expression to sum over a set

  • Use the sort operator

Problem¶

../_images/stochastic_packing.svg

The stochastic packing problem is defined as a set of items that need to be grouped into bins. Each bin can contain any item, but the items have random weights. Here, randomness is represented by scenarios, where one scenario defines the weights of all items. Given a distribution of items in bins, each scenario yields a bin of maximum weight.

Now, which objective is most appropriate to minimize this stochastic maximum weight?

Minimizing the average can hide risky scenarios while minimizing the worst-case might be too pessimistic. A usual compromise to build a robust scenario is to optimize on a given percentile. It is what is done, minimizing the 90 th percentile of makespans. Thanks to the sort operator, such a nonlinear criterion is straightforward to model using Hexaly Optimizer.

Download the example


Data¶

For this example, we generate instances at runtime: first a uniform distribution is picked for each item, then for each scenario, the weight of each item is independently sampled from the corresponding uniform distribution.

Program¶

We use a set decision variable to represent the set of items contained in a bin. Those sets are constrained to form a partition as each item must be present in exactly one bin.

We then compute and store in the scenarioMaxWeight array the maximum weight corresponding to each scenario. To do so, we first need to compute the total weight of each bin as a sum over a set and therefore need to define a lambda expression.

We can then sort the array of maximum weights and access our objective function: its 9 th decile.

Execution:
hexaly stochastic_packing.hxm [hxTimeLimit] [solFileName=solution.txt]
use random;

/* Generate instance data */
function input() {
    nbItems = 10;
    nbBins = 2;
    nbScenarios = 3;
    rngSeed = 42;

    // Pick random parameters for each item distribution
    rng = random.create(rngSeed);
    itemsMin[i in 0...nbItems] = rng.next(10, 101);
    itemsMax[i in 0...nbItems] = itemsMin[i] + rng.next(0, 51);

    // Sample the distributions to generate the scenarios
    scenarioItemWeights[i in 0...nbScenarios][j in 0...nbItems] = 
            rng.next(itemsMin[j], itemsMax[j] + 1);
}

/* Declare the optimization model */
function model() {
    // Set decisions: bins[k] represents the items in bin k
    bins[k in 0...nbBins] <- set(nbItems);

    // Each item must be in one bin and one bin only
    constraint partition[k in 0...nbBins](bins[k]);

    // Compute max weight for each scenario
    scenarioMaxWeight[m in 0...nbScenarios] <- max[k in 0...nbBins](
            sum(bins[k], i => scenarioItemWeights[m][i]));

    // Compute the 9th decile of scenario max weights
    stochasticMaxWeight <- sort(scenarioMaxWeight)[ceil(0.9 * (nbScenarios - 1))];

    minimize stochasticMaxWeight;
}

// Parametrize the solver
function param() {
    if (hxTimeLimit == nil) hxTimeLimit = 2;
}

/* Write the solution */
function output() {
    println();
    println("Scenario item weights:");
    for [i in 0...nbScenarios] {
        print(i + ": [");
        for [j in 0...nbItems]
            print(scenarioItemWeights[i][j] + (j == nbItems - 1 ? "" : ", "));
        println("]");
    }

    println();
    println("Bins:");
    for [k in 0...nbBins]
        println(k + ": " + bins[k].value);
}
Execution (Windows)
set PYTHONPATH=%HX_HOME%\bin\python
python stochastic_packing.py
Execution (Linux)
export PYTHONPATH=/opt/hexaly_13_0/bin/python
python stochastic_packing.py
from __future__ import print_function
import random
import math
import hexaly.optimizer


def generate_scenarios(nb_items, nb_scenarios, rng_seed):
    random.seed(rng_seed)

    # Pick random parameters for each item distribution
    items_dist = []
    for _ in range(nb_items):
        item_min = random.randint(10, 100)
        item_max = item_min + random.randint(0, 50)
        items_dist.append((item_min, item_max))

    # Sample the distributions to generate the scenarios
    scenario_item_weights = [[random.randint(*dist) for dist in items_dist]
                               for _ in range(nb_scenarios)]
    return scenario_item_weights


def main(nb_items, nb_bins, nb_scenarios, seed, time_limit):
    # Generate instance data
    scenario_item_weights_data = generate_scenarios(nb_items, nb_scenarios, seed)

    with hexaly.optimizer.HexalyOptimizer() as optimizer:
        #
        # Declare the optimization model
        #
        model = optimizer.model

        # Set decisions: bins[k] represents the items in bin k
        bins = [model.set(nb_items) for _ in range(nb_bins)]

        # Each item must be in one bin and one bin only
        model.constraint(model.partition(bins))

        scenarios_item_weights = model.array(scenario_item_weights_data)

        # Compute max weight for each scenario
        scenarios_max_weights = model.array(
            model.max(
                model.sum(bin,
                          model.lambda_function(
                              lambda i:
                                  model.at(scenarios_item_weights, k, i)))
                for bin in bins) for k in range(nb_scenarios))

        # Compute the 9th decile of scenario max weights
        stochastic_max_weight = \
            model.sort(scenarios_max_weights)[int(math.ceil(0.9 * (nb_scenarios - 1)))]

        model.minimize(stochastic_max_weight)
        model.close()

        # Parameterize the optimizer
        optimizer.param.time_limit = time_limit

        optimizer.solve()

        #
        # Write the solution
        #
        print()
        print("Scenario item weights:")
        for i, scenario in enumerate(scenario_item_weights_data):
            print(i, ': ', scenario, sep='')

        print()
        print("Bins:")
        for k, bin in enumerate(bins):
            print(k, ': ', bin.value, sep='')


if __name__ == '__main__':
    nb_items = 10
    nb_bins = 2
    nb_scenarios = 3
    rng_seed = 42
    time_limit = 2

    main(
        nb_items,
        nb_bins,
        nb_scenarios,
        rng_seed,
        time_limit
    )
Compilation / Execution (Windows)
cl /EHsc stochastic_packing.cpp -I%HX_HOME%\include /link %HX_HOME%\bin\hexaly130.lib
stochastic_packing
Compilation / Execution (Linux)
g++ stochastic_packing.cpp -I/opt/hexaly_13_0/include -lhexaly130 -lpthread -o stochastic_packing
stochastic_packing
#include "optimizer/hexalyoptimizer.h"
#include <cmath>
#include <iostream>
#include <random>
#include <vector>

using namespace hexaly;

class StochasticPacking {
private:
    // Number of items
    int nbItems;
    // Number of bins
    int nbBins;
    // Number of scenarios
    int nbScenarios;
    // For each scenario, the weight of each item
    std::vector<std::vector<int>> scenarioItemWeights;

    // Hexaly Optimizer
    HexalyOptimizer optimizer;
    // Decision variable for the assignment of items
    std::vector<HxExpression> bins;
    // For each scenario, the corresponding maximum weight
    std::vector<HxExpression> scenarioMaxWeight;
    // Objective = minimize the 9th decile of all possible max weights
    HxExpression stochasticMaxWeight;

    void generateScenarios(unsigned int rngSeed) {
        std::mt19937 rng(rngSeed);
        std::uniform_int_distribution<int> distMin(10, 100);
        std::uniform_int_distribution<int> distDelta(0, 50);

        // Pick random parameters for each item distribution
        std::vector<std::uniform_int_distribution<int>> itemsDists;
        for (int i = 0; i < nbItems; ++i) {
            int min = distMin(rng);
            int max = min + distDelta(rng);
            itemsDists.emplace_back(min, max);
        }

        // Sample the distributions to generate the scenarios
        for (int i = 0; i < nbScenarios; ++i) {
            for (int j = 0; j < nbItems; ++j) {
                scenarioItemWeights[i][j] = itemsDists[j](rng);
            }
        }
    }

public:
    StochasticPacking(int nbItems, int nbBins, int nbScenarios, unsigned int seed)
        : nbItems(nbItems), nbBins(nbBins), nbScenarios(nbScenarios),
          scenarioItemWeights(nbScenarios, std::vector<int>(nbItems)), optimizer() {
        generateScenarios(seed);
    }

    void solve(int timeLimit) {
        // Declare the optimization model
        HxModel model = optimizer.getModel();

        bins.resize(nbBins);
        scenarioMaxWeight.resize(nbScenarios);

        // Set decisions: bins[k] represents the items in bin k
        for (int k = 0; k < nbBins; ++k) {
            bins[k] = model.setVar(nbItems);
        }

        // Each item must be in one bin and one bin only
        model.constraint(model.partition(bins.begin(), bins.end()));

        // Compute max weight for each scenario
        for (int m = 0; m < nbScenarios; ++m) {
            HxExpression scenario = model.array(scenarioItemWeights[m].begin(), scenarioItemWeights[m].end());
            HxExpression weightLambda = model.createLambdaFunction([&](HxExpression i) { return scenario[i]; });
            std::vector<HxExpression> binWeights(nbBins);

            for (int k = 0; k < nbBins; ++k) {
                binWeights[k] = model.sum(bins[k], weightLambda);
            }
            scenarioMaxWeight[m] = model.max(binWeights.begin(), binWeights.end());
        }

        // Compute the 9th decile of scenario max weights
        HxExpression scenarioMaxWeightArray = model.array(scenarioMaxWeight.begin(), scenarioMaxWeight.end());
        HxExpression sortedScenarioMaxWeight = model.sort(scenarioMaxWeightArray);
        stochasticMaxWeight = sortedScenarioMaxWeight[(int)std::ceil(0.9 * (nbScenarios - 1))];

        model.minimize(stochasticMaxWeight);
        model.close();

        // Parametrize the optimizer
        optimizer.getParam().setTimeLimit(timeLimit);

        optimizer.solve();
    }

    /* Write the solution */
    void writeSolution(std::ostream& os) const {
        os << "\nScenario item weights:\n";
        for (int i = 0; i < nbScenarios; ++i) {
            os << i << ": [";
            for (int j = 0; j < scenarioItemWeights[i].size(); ++j) {
                os << scenarioItemWeights[i][j] << (j == scenarioItemWeights[i].size() - 1 ? "" : ", ");
            }
            os << "]\n";
        }

        os << "\nBins:\n";
        for (int m = 0; m < nbBins; ++m) {
            os << m << ": { ";
            HxCollection items = bins[m].getCollectionValue();
            for (int i = 0; i < items.count(); ++i) {
                os << items[i] << (i == items.count() - 1 ? " " : ", ");
            }
            os << "}\n";
        }
    }
};

int main(int argc, char** argv) {
    int nbItems = 10;
    int nbBins = 2;
    int nbScenarios = 3;
    int rngSeed = 42;
    int timeLimit = 2;

    try {
        StochasticPacking model(nbItems, nbBins, nbScenarios, rngSeed);
        model.solve(timeLimit);
        model.writeSolution(std::cout);
        return 0;
    } catch (const std::exception& e) {
        std::cerr << "An error occurred: " << e.what() << std::endl;
        return 1;
    }
}
Compilation/Execution (Windows)
copy %HX_HOME%\bin\Hexaly.NET.dll .
csc StochasticPacking.cs /reference:Hexaly.NET.dll
StochasticPacking
using System;

using Hexaly.Optimizer;

public class StochasticPacking : IDisposable
{
    // Number of items
    int nbItems;

    // Number of bins
    int nbBins;

    // Number of scenarios
    int nbScenarios;

    // For each scenario, the weight of each item
    int[][] scenarioItemWeights;

    // Hexaly Optimizer
    HexalyOptimizer optimizer;

    // Decision variable for the assignment of items
    HxExpression[] bins;

    // For each scenario, the corresponding maximum weight
    HxExpression[] scenarioMaxWeight;

    // Objective = minimize the 9th decile of all possible max weights
    HxExpression stochasticMaxWeight;

    private void generateScenarios(int rngSeed)
    {
        Random rng = new Random(rngSeed);

        // Pick random parameters for each item distribution
        int[] itemsMin = new int[nbItems];
        int[] itemsMax = new int[nbItems];
        for (int i = 0; i < nbItems; ++i)
        {
            itemsMin[i] = rng.Next(10, 101);
            itemsMax[i] = itemsMin[i] + rng.Next(51);
        }

        // Sample the distributions to generate the scenarios
        scenarioItemWeights = new int[nbScenarios][];
        for (int i = 0; i < nbScenarios; ++i)
        {
            scenarioItemWeights[i] = new int[nbItems];
            for (int j = 0; j < nbItems; ++j)
                scenarioItemWeights[i][j] = rng.Next(itemsMin[i], itemsMax[i] + 1);
        }
    }

    public StochasticPacking(int nbItems, int nbBins, int nbScenarios, int rngSeed)
    {
        optimizer = new HexalyOptimizer();
        this.nbItems = nbItems;
        this.nbBins = nbBins;
        this.nbScenarios = nbScenarios;
        generateScenarios(rngSeed);
    }

    public void Dispose()
    {
        if (optimizer != null)
            optimizer.Dispose();
    }

    void Solve(int limit)
    {
        // Declare the optimization model
        HxModel model = optimizer.GetModel();

        bins = new HxExpression[nbBins];
        scenarioMaxWeight = new HxExpression[nbScenarios];

        // Set decisions: bins[k] represents the items in bin k
        for (int k = 0; k < nbBins; ++k)
            bins[k] = model.Set(nbItems);

        // Each item must be in one bin and one bin only
        model.Constraint(model.Partition(bins));

        // Compute max weight for each scenario
        for (int m = 0; m < nbScenarios; ++m)
        {
            HxExpression scenario = model.Array(scenarioItemWeights[m]);
            HxExpression weightLambda = model.LambdaFunction(i => scenario[i]);
            HxExpression[] binWeights = new HxExpression[nbBins];

            for (int k = 0; k < nbBins; ++k)
                binWeights[k] = model.Sum(bins[k], weightLambda);
            scenarioMaxWeight[m] = model.Max(binWeights);
        }

        // Compute the 9th decile of scenario max weights
        HxExpression scenarioMaxWeightArray = model.Array(scenarioMaxWeight);
        HxExpression sortedScenarioMaxWeight = model.Sort(scenarioMaxWeightArray);
        stochasticMaxWeight = sortedScenarioMaxWeight[(int)Math.Ceiling(0.9 * (nbScenarios - 1))];

        model.Minimize(stochasticMaxWeight);
        model.Close();

        // Parametrize the optimizer
        optimizer.GetParam().SetTimeLimit(limit);

        optimizer.Solve();
    }

    /* Write the solution */
    private void WriteSolution()
    {
        Console.WriteLine();
        Console.WriteLine("Scenario item weights:");
        for (int i = 0; i < nbScenarios; ++i)
        {
            Console.Write(i + ": [");
            for (int j = 0; j < nbItems; ++j)
                Console.Write(scenarioItemWeights[i][j] + (j == nbItems - 1 ? "" : ", "));
            Console.WriteLine("]");
        }

        Console.WriteLine();
        Console.WriteLine("Bins:");

        for (int m = 0; m < nbBins; ++m)
        {
            Console.Write(m + ": { ");
            HxCollection items = bins[m].GetCollectionValue();
            for (int i = 0; i < items.Count(); ++i)
                Console.Write(items.Get(i) + (i == items.Count() - 1 ? " " : ", "));
            Console.WriteLine("}");
        }
    }

    public static void Main(string[] args)
    {
        int nbItems = 10;
        int nbBins = 2;
        int nbScenarios = 3;
        int rngSeed = 43;
        int timeLimit = 2;

        using (
            StochasticPacking model = new StochasticPacking(
                nbItems,
                nbBins,
                nbScenarios,
                rngSeed
            )
        )
        {
            model.Solve(timeLimit);
            model.WriteSolution();
        }
    }
}
Compilation / Execution (Windows)
javac StochasticPacking.java -cp %HX_HOME%\bin\hexaly.jar
java -cp %HX_HOME%\bin\hexaly.jar;. StochasticPacking
Compilation/Execution (Linux)
javac StochasticPacking.java -cp /opt/hexaly_13_0/bin/hexaly.jar
java -cp /opt/hexaly_13_0/bin/hexaly.jar:. StochasticPacking
import java.util.Random;

import com.hexaly.optimizer.*;

public class StochasticPacking {
    // Number of items
    private int nbItems;

    // Number of bins
    private int nbBins;

    // Number of scenarios
    private int nbScenarios;

    // For each scenario, the weight of each item
    private int[][] scenarioItemWeights;

    // Hexaly Optimizer
    private final HexalyOptimizer optimizer;

    // Decision variable for the assignment of items
    private HxExpression[] bins;

    // For each scenario, the corresponding max weight
    private HxExpression[] scenarioMaxWeight;

    // Objective = minimize the 9th decile of all possible max weights
    private HxExpression stochasticMaxWeight;

    private void generateScenarios(int rngSeed) {
        Random rng = new Random(rngSeed);

        // Pick random parameters for each item distribution
        int[] itemsMin = new int[nbItems];
        int[] itemsMax = new int[nbItems];
        for (int i = 0; i < nbItems; ++i) {
            itemsMin[i] = 10 + rng.nextInt(91);
            itemsMax[i] = itemsMin[i] + rng.nextInt(51);
        }

        // Sample the distributions to generate the scenarios
        scenarioItemWeights = new int[nbScenarios][nbItems];
        for (int i = 0; i < nbScenarios; ++i) {
            for (int j = 0; j < nbItems; ++j) {
                scenarioItemWeights[i][j] = itemsMin[j] + rng.nextInt(itemsMax[i] - itemsMin[i] + 1);
            }
        }
    }

    private StochasticPacking(HexalyOptimizer optimizer, int nbItems, int nbBins, int nbScenarios, int rngSeed) {
        this.optimizer = optimizer;
        this.nbItems = nbItems;
        this.nbBins = nbBins;
        this.nbScenarios = nbScenarios;
        generateScenarios(rngSeed);
    }

    private void solve(int limit) {
        // Declare the optimization model
        HxModel model = optimizer.getModel();

        bins = new HxExpression[nbBins];
        scenarioMaxWeight = new HxExpression[nbScenarios];

        // Set decisions: bins[k] represents the items in bin k
        for (int k = 0; k < nbBins; ++k) {
            bins[k] = model.setVar(nbItems);
        }

        // Each item must be in one bin and one bin only
        model.constraint(model.partition(bins));

        // Compute max weight for each scenario
        for (int m = 0; m < nbScenarios; ++m) {
            HxExpression scenario = model.array(scenarioItemWeights[m]);
            HxExpression weightLambda = model.lambdaFunction(i -> model.at(scenario, i));
            HxExpression[] binWeights = new HxExpression[nbBins];

            for (int k = 0; k < nbBins; ++k) {
                binWeights[k] = model.sum(bins[k], weightLambda);
            }
            scenarioMaxWeight[m] = model.max(binWeights);
        }

        // Compute the 9th decile of scenario makespans
        HxExpression scenarioMaxWeightArray = model.array(scenarioMaxWeight);
        HxExpression sortedScenarioMaxWeight = model.sort(scenarioMaxWeightArray);
        stochasticMaxWeight = model.at(sortedScenarioMaxWeight, (int) Math.ceil(0.9 * (nbScenarios - 1)));

        model.minimize(stochasticMaxWeight);
        model.close();

        // Parametrize the optimizer
        optimizer.getParam().setTimeLimit(limit);

        optimizer.solve();
    }

    /* Write the solution */
    private void writeSolution() {
        System.out.println();
        System.out.println("Scenario item weights:");
        for (int i = 0; i < nbScenarios; ++i) {
            System.out.print("" + i + ": [");
            for (int j = 0; j < nbItems; ++j) {
                System.out.print("" + scenarioItemWeights[i][j] + (j == nbItems - 1 ? "" : ", "));
            }
            System.out.println("]");
        }

        System.out.println();
        System.out.println("Bins:");

        for (int m = 0; m < nbBins; ++m) {
            System.out.print("" + m + ": { ");
            HxCollection items = bins[m].getCollectionValue();
            for (int i = 0; i < items.count(); ++i) {
                System.out.print("" + items.get(i) + (i == items.count() - 1 ? " " : ", "));
            }
            System.out.println("}");
        }
    }

    public static void main(String[] args) {
        try (HexalyOptimizer optimizer = new HexalyOptimizer()) {
            int nbItems = 10;
            int nbBins = 2;
            int nbScenarios = 3;
            int rngSeed = 42;
            int timeLimit = 2;

            StochasticPacking model = new StochasticPacking(optimizer, nbItems, nbBins, nbScenarios,
                rngSeed);

            model.solve(timeLimit);
            model.writeSolution();
        } catch (Exception ex) {
            System.err.println(ex);
            ex.printStackTrace();
            System.exit(1);
        }
    }
};