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
- Use set decision variables to model the contents of the bins
- Define a lambda function to compute the total weight of a bin
- Use the ‘sort’ operator to compute the 90th percentile
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\pythonpython stochastic_packing.py
- Execution (Linux)
-
export PYTHONPATH=/opt/hexaly_13_0/bin/pythonpython 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.libstochastic_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.dllStochasticPacking
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.jarjava -cp %HX_HOME%\bin\hexaly.jar;. StochasticPacking
- Compilation / Execution (Linux)
-
javac StochasticPacking.java -cp /opt/hexaly_13_0/bin/hexaly.jarjava -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);
}
}
};