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:
localsolver stochastic_packing.lsp [lsTimeLimit=n] [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 (lsTimeLimit == nil) lsTimeLimit = 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=%LS_HOME%\bin\python
python stochastic_packing.py
Execution (Linux)
export PYTHONPATH=/opt/localsolver_12_5/bin/python
python stochastic_packing.py
from __future__ import print_function
import random
import math
import localsolver


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 localsolver.LocalSolver() as ls:
        #
        # Declare the optimization model
        #
        model = ls.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 solver
        ls.param.time_limit = time_limit

        ls.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%LS_HOME%\include /link %LS_HOME%\bin\localsolver125.lib
stochastic_packing
Compilation / Execution (Linux)
g++ stochastic_packing.cpp -I/opt/localsolver_12_5/include -llocalsolver125 -lpthread -o stochastic_packing
stochastic_packing
#include "localsolver.h"
#include <cmath>
#include <iostream>
#include <random>
#include <vector>

using namespace localsolver;

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;

    // LocalSolver
    LocalSolver localsolver;
    // Decision variable for the assignment of items
    std::vector<LSExpression> bins;
    // For each scenario, the corresponding maximum weight
    std::vector<LSExpression> scenarioMaxWeight;
    // Objective = minimize the 9th decile of all possible max weights
    LSExpression 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)), localsolver() {
        generateScenarios(seed);
    }

    void solve(int timeLimit) {
        // Declare the optimization model
        LSModel model = localsolver.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) {
            LSExpression scenario = model.array(scenarioItemWeights[m].begin(), scenarioItemWeights[m].end());
            LSExpression weightLambda = model.createLambdaFunction([&](LSExpression i) { return scenario[i]; });
            std::vector<LSExpression> 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
        LSExpression scenarioMaxWeightArray = model.array(scenarioMaxWeight.begin(), scenarioMaxWeight.end());
        LSExpression sortedScenarioMaxWeight = model.sort(scenarioMaxWeightArray);
        stochasticMaxWeight = sortedScenarioMaxWeight[(int)std::ceil(0.9 * (nbScenarios - 1))];

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

        // Parametrize the solver
        localsolver.getParam().setTimeLimit(timeLimit);

        localsolver.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 << ": { ";
            LSCollection 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 %LS_HOME%\bin\localsolvernet.dll .
csc StochasticPacking.cs /reference:localsolvernet.dll
StochasticPacking
using System;

using localsolver;

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;

    // LocalSolver
    LocalSolver localsolver;

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

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

    // Objective = minimize the 9th decile of all possible max weights
    LSExpression 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)
    {
        localsolver = new LocalSolver();
        this.nbItems = nbItems;
        this.nbBins = nbBins;
        this.nbScenarios = nbScenarios;
        generateScenarios(rngSeed);
    }

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

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

        bins = new LSExpression[nbBins];
        scenarioMaxWeight = new LSExpression[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)
        {
            LSExpression scenario = model.Array(scenarioItemWeights[m]);
            LSExpression weightLambda = model.LambdaFunction(i => scenario[i]);
            LSExpression[] binWeights = new LSExpression[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
        LSExpression scenarioMaxWeightArray = model.Array(scenarioMaxWeight);
        LSExpression sortedScenarioMaxWeight = model.Sort(scenarioMaxWeightArray);
        stochasticMaxWeight = sortedScenarioMaxWeight[(int)Math.Ceiling(0.9 * (nbScenarios - 1))];

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

        // Parametrize the solver
        localsolver.GetParam().SetTimeLimit(limit);

        localsolver.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 + ": { ");
            LSCollection 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 %LS_HOME%\bin\localsolver.jar
java -cp %LS_HOME%\bin\localsolver.jar;. StochasticPacking
Compilation/Execution (Linux)
javac StochasticPacking.java -cp /opt/localsolver_12_5/bin/localsolver.jar
java -cp /opt/localsolver_12_5/bin/localsolver.jar:. StochasticPacking
import java.util.Random;

import localsolver.*;

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;

    // LocalSolver
    private final LocalSolver localsolver;

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

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

    // Objective = minimize the 9th decile of all possible max weights
    private LSExpression 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(LocalSolver localsolver, int nbItems, int nbBins, int nbScenarios, int rngSeed) {
        this.localsolver = localsolver;
        this.nbItems = nbItems;
        this.nbBins = nbBins;
        this.nbScenarios = nbScenarios;
        generateScenarios(rngSeed);
    }

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

        bins = new LSExpression[nbBins];
        scenarioMaxWeight = new LSExpression[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) {
            LSExpression scenario = model.array(scenarioItemWeights[m]);
            LSExpression weightLambda = model.lambdaFunction(i -> model.at(scenario, i));
            LSExpression[] binWeights = new LSExpression[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
        LSExpression scenarioMaxWeightArray = model.array(scenarioMaxWeight);
        LSExpression sortedScenarioMaxWeight = model.sort(scenarioMaxWeightArray);
        stochasticMaxWeight = model.at(sortedScenarioMaxWeight, (int) Math.ceil(0.9 * (nbScenarios - 1)));

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

        // Parametrize the solver
        localsolver.getParam().setTimeLimit(limit);

        localsolver.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 + ": { ");
            LSCollection 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 (LocalSolver localsolver = new LocalSolver()) {
            int nbItems = 10;
            int nbBins = 2;
            int nbScenarios = 3;
            int rngSeed = 42;
            int timeLimit = 2;

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

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