Stochastic Packing Problem

Problem

In the Stochastic Packing Problem, a set of items must be grouped into bins. There are no capacity constraints, and the bins can contain any item. The stochastic aspect of the problem comes from the fact that the items have random weights. In this example, randomness is represented by different scenarios, corresponding to different possible item weights. For a given distribution of items into bins, each scenario is associated with a bin of maximum weight.

Determining the most relevant objective function to minimize this stochastic maximum weight can be a challenge. Indeed, minimizing the average maximum weight can hide risky scenarios, while minimizing the maximum weight in the worst-case scenario might be too pessimistic. A usual compromise to build a robust scenario is to optimize on a given percentile. Consequently, the objective function we use in this example is the minimization of the 90th percentile of maximum weights among all scenarios.

Principles learned

Data

For this example, we generate random instances at runtime. We start by picking a uniform distribution for each item. Then, for each scenario, the weight of each item is independently sampled from the corresponding uniform distribution.

Program

The Hexaly model for the Stochastic Packing Problem uses set decision variables to represent the set of items contained in each bin. To ensure that each item is assigned to exactly one bin, we constrain the sets to form a partition.

We compute the total weight of a bin in each scenario using a variadic ‘sum’ operator on the set and a lambda function returning the weight associated with any item index in this scenario. Note that the number of terms in this sum varies during the search, along with the size of the set. We can then compute the maximum bin weight associated with each scenario.

Finally, we use the ‘sort’ operator to sort the array of maximum weights in ascending order. We can then easily access the value of the objective function, which is the 9th decile in the sorted array.

Execution
hexaly stochastic_packing.hxm [hxTimeLimit=] [solFileName=]
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);
        }
    }
};