Stochastic Scheduling

Principles learned

  • Use list decision variables

  • Use a lambda expression to sum over a set

  • Use the sort operator

Problem

../_images/stochastic_scheduling.svg

The Stochastic Scheduling Problem is defined as a set of tasks that has to be processed by machines. Any machine can independently process all tasks, but tasks have random processing times. Here, the randomness is represented by scenarios, where a scenario defines the processing times of all the tasks. Given a processing schedule, each scenario yields a makespan: the time when all tasks have been processed.

Now, which objective is most appropriate to minimize the stochastic makespan?

Minimizing the average makespan 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 LocalSolver.

Download the example


Data

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

Program

We use a set decision variable to represent the set of tasks processed by each machine. Those sets are constrained to form a partition as each task must be processed by exactly one machine.

We then compute and store in the scenarioMakespan array the makespan corresponding to each scenario. To do so, we first need to compute the total processing time of each machine as a sum over a set and therefore need to define a lambda expression.

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

Execution:
localsolver stochastic_scheduling.lsp [lsTimeLimit=n] [solFileName=solution.txt]
/********** stochastic_scheduling.lsp **********/

use random;

// Generate instance data.
function input() {
    nbTasks = 10;
    nbMachines = 2;
    nbScenarios = 3;
    rngSeed = 42;

    // Pick random parameters for each task distribution
    rng = random.create(rngSeed);
    taskMin[0..nbTasks - 1] = rng.next(10, 101);
    taskMax[i in 0..nbTasks - 1] = taskMin[i] + rng.next(0, 51);

    // Sample the distributions to generate the scenarios
    scenarioTaskDurations[0..nbScenarios - 1][j in 0..nbTasks - 1] = rng.next(taskMin[j], taskMax[j] + 1);
}

// Declares the optimization model.
function model() {
    // Set decisions: machineTasks[k] represents the tasks processed my machine k
    machineTasks[k in 0..nbMachines - 1] <- set(nbTasks);

    // Each task must be processed by exactly one machine
    constraint partition[k in 0..nbMachines - 1](machineTasks[k]);

    // Compute makespan of each scenario
    scenarioMakespan[m in 0..nbScenarios - 1] <-
        max[k in 0..nbMachines - 1] (
            sum(machineTasks[k],
                i => scenarioTaskDurations[m][i]
            )
        );

    // Compute the 90th percentile of scenario makespans.
    stochasticMakespan <- sort(scenarioMakespan)[ceil(0.9 * (nbScenarios - 1))];

    minimize stochasticMakespan;
}

// Parameterizes the solver.
function param() {
    if (lsTimeLimit == nil) lsTimeLimit = 2;
}

function output() {
    println();
    println("Scenario task durations:");
    for [i in 0..nbScenarios - 1] {
        print(i + ": [");
        for [j in 0..nbTasks - 1]
            print(scenarioTaskDurations[i][j] + (j == nbTasks - 1 ? "" : ", "));
        println("]");
    }

    println();
    println("Machines:");
    for [k in 0..nbMachines - 1]
        println(k + ": " + machineTasks[k].value);
}
Execution (Windows)
set PYTHONPATH=%LS_HOME%\bin\python27\
python stochastic_scheduling.py
Execution (Linux)
export PYTHONPATH=/opt/localsolver_XXX/bin/python27/
python stochastic_scheduling.py
########## stochastic_scheduling.py ##########

from __future__ import print_function
import random
import math

import localsolver

def generate_scenarios(nb_tasks, nb_scenarios, rng_seed):
    random.seed(rng_seed)
    
    # Pick random parameters for each task distribution
    tasks_dist = []
    for _ in range(nb_tasks):
        task_min = random.randint(10, 100)
        task_max = task_min + random.randint(0, 50)
        tasks_dist.append((task_min, task_max))
    
    # Sample the distributions to generate the scenarios
    scenario_task_durations = [
        [
            random.randint(*dist) for dist in tasks_dist
        ] for _ in range(nb_scenarios)
    ]
    return scenario_task_durations

def main(nb_tasks, nb_machines, nb_scenarios, seed, time_limit):
    # Generate instance data
    scenario_task_durations = generate_scenarios(nb_tasks, nb_scenarios, seed)

    with localsolver.LocalSolver() as ls:
        # Declares the optimization model
        model = ls.model

        # Set decisions: machineTasks[k] represents the tasks processed my machine k
        machine_tasks = [model.set(nb_tasks) for _ in range(nb_machines)]

        # Each task must be processed by exactly one machine
        model.constraint(model.partition(machine_tasks))

        scenarios_task_durations_array = model.array(scenario_task_durations)
        
        # Compute makespan of each scenario
        scenarios_makespan = model.array(
            model.max(
                model.sum(
                    machine,
                    model.lambda_function(
                        lambda i : model.at(scenarios_task_durations_array, k, i)
                    )
                )
                for machine in machine_tasks
            )
            for k in range(nb_scenarios)
        )

        # Compute the 90th percentile of scenario makespans.
        stochastic_makespan = \
            model.sort(scenarios_makespan)[int(math.ceil(0.9 * (nb_scenarios - 1)))]

        model.minimize(stochastic_makespan)
        model.close()

        # Parameterizes the solver
        ls.param.time_limit = time_limit

        ls.solve()

        print()
        print("Scenario task durations:")
        for i, scenario in enumerate(scenario_task_durations):
            print(i, ': ', scenario, sep='')

        print()
        print("Machines:")
        for k, machine in enumerate(machine_tasks):
            print(k, ': ', machine.value, sep='')


if __name__ == '__main__':
    nb_tasks = 10
    nb_machines = 2
    nb_scenarios = 3
    rng_seed = 42
    time_limit = 2

    main(
        nb_tasks,
        nb_machines,
        nb_scenarios,
        rng_seed,
        time_limit
    )
Compilation / Execution (Windows)
cl /EHsc stochastic_scheduling.cpp -I%LS_HOME%\include /link %LS_HOME%\bin\localsolver.dll.lib
stochastic_scheduling
Compilation / Execution (Linux)
g++ stochastic_scheduling.cpp -I/opt/localsolver_XXX/include -llocalsolver -lpthread -o stochastic_scheduling
stochastic_scheduling
/********** stochastic_scheduling.cpp **********/

#include <iostream>
#include <random>
#include <vector>
#include <cmath>

#include "localsolver.h"

using namespace localsolver;

class StochasticScheduling {
private:
    //Number of tasks
    int nbTasks;
    //Number of machines
    int nbMachines;
    //Number of scenarios
    int nbScenarios;
    //For each scenario, the processing time of each task
    std::vector<std::vector<int> > scenarioTaskDurations;

    // LocalSolver
    LocalSolver localsolver;
    // Decision variable for the assignment of tasks
    std::vector<LSExpression> machineTasks;
    // For each scenario, the corresponding makespan
    std::vector<LSExpression> scenarioMakespan;
    // Objective = minimize the 9th decile of all possible makespans
    LSExpression stochasticMakespan;

    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 task distribution
        std::vector<std::uniform_int_distribution<int> > tasksDists;
        for (int i = 0; i < nbTasks; i++) {
            int min = distMin(rng);
            int max = min + distDelta(rng);
            tasksDists.emplace_back(min, max);
        }

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

public:

    StochasticScheduling(int nbTasks, int nbMachines, int nbScenarios, unsigned int seed)
        : nbTasks(nbTasks),
          nbMachines(nbMachines),
          nbScenarios(nbScenarios),
          scenarioTaskDurations(nbScenarios, std::vector<int>(nbTasks)),
          localsolver() {
        generateScenarios(seed);
    }

    void solve(int timeLimit) {
        // Declares the optimization model.
        LSModel model = localsolver.getModel();

        machineTasks.resize(nbMachines);
        scenarioMakespan.resize(nbScenarios);

        // Set decisions: machineTasks[k] represents the tasks processed my machine k
        for (int k = 0; k < nbMachines; k++) {
            machineTasks[k] = model.setVar(nbTasks);
        }

        // Each task must be processed by exactly one machine
        model.constraint(model.partition(machineTasks.begin(), machineTasks.end()));

        // Compute makespan of each scenario
        for (int m = 0; m < nbScenarios; m++) {
            LSExpression scenario = model.array(scenarioTaskDurations[m].begin(), scenarioTaskDurations[m].end());
            LSExpression taskSelector = model.createLambdaFunction([&](LSExpression i){ return scenario[i]; });
            std::vector<LSExpression> machinesMakespan(nbMachines);
            
            for (int k = 0; k < nbMachines; k++) {
                machinesMakespan[k] = model.sum(
                    machineTasks[k], 
                    taskSelector
                );
            }
            scenarioMakespan[m] = model.max(machinesMakespan.begin(), machinesMakespan.end());
        }

        // Compute the 90th percentile of scenario makespans.
        LSExpression scenarioMakespanArray = model.array(scenarioMakespan.begin(), scenarioMakespan.end());
        LSExpression sortedScenarioMakespan = model.sort(scenarioMakespanArray);
        stochasticMakespan = sortedScenarioMakespan[(int)std::ceil(0.9 * (nbScenarios - 1))];

        model.minimize(stochasticMakespan);
        model.close();

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

        localsolver.solve();
    }

    void writeSolution(std::ostream &os) const {
        os << "\nScenario task durations:\n";
        for (int i = 0; i < nbScenarios; i++) {
            os << i << ": [";
            for (int j = 0; j < scenarioTaskDurations[i].size(); j++) {
                os << scenarioTaskDurations[i][j] << (j == scenarioTaskDurations[i].size() - 1 ? "" : ", ");
            }
            os << "]\n";
        }


        os << "\nMachines:\n";
        for (int m = 0; m < nbMachines; m++) {
            os << m << ": { ";
            LSCollection tasks = machineTasks[m].getCollectionValue();
            for (int i = 0; i < tasks.count(); i++) {
                os << tasks[i] << (i == tasks.count() - 1 ? " " : ", ");
            }
            os << "}\n";
        }
    }
};

int main(int argc, char** argv) {
    int nbTasks = 10;
    int nbMachines = 2;
    int nbScenarios = 3;
    int rngSeed = 42;
    int timeLimit = 2;
    
    try {
        StochasticScheduling model(
            nbTasks,
            nbMachines,
            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)
javac StochasticScheduling.java -cp %LS_HOME%\bin\localsolver.jar
java -cp %LS_HOME%\bin\localsolver.jar;. StochasticScheduling
Compilation/Execution (Linux)
javac StochasticScheduling.java -cp /opt/localsolver_XXX/bin/localsolver.jar
java -cp /opt/localsolver_XXX/bin/localsolver.jar:. StochasticScheduling
/********** stochastic_scheduling.java **********/

import java.io.*;
import java.util.Random;

import localsolver.*;

public class StochasticScheduling {
    //Number of tasks
    private int nbTasks;
   
    //Number of machines
    private int nbMachines;
   
    //Number of scenarios
    private int nbScenarios;

    //For each scenario, the processing time of each task
    private int[][] scenarioTaskDurations;

    // Solver
    private final LocalSolver localsolver;

    // Decision variable for the assignment of tasks
    private LSExpression[] machineTasks;
    
    // For each scenario, the corresponding makespan
    private LSExpression[] scenarioMakespan;
    
    // Objective = minimize the 9th decile of all possible makespans
    private LSExpression stochasticMakespan;

    private void generateScenarios(int rngSeed) {
        Random rng = new Random(rngSeed);
        
        // Pick random parameters for each task distribution
        int[] tasksMin = new int[nbTasks];
        int[] tasksMax = new int[nbTasks];
        for (int i = 0; i < nbTasks; i++) {
            tasksMin[i] = 10 + rng.nextInt(91);
            tasksMax[i] = tasksMin[i] + rng.nextInt(51);
        }

        // Sample the distributions to generate the scenarios
        scenarioTaskDurations = new int[nbScenarios][nbTasks];
        for (int i = 0; i < nbScenarios; i++) {
            for (int j = 0; j < nbTasks; ++j) {
                scenarioTaskDurations[i][j] = tasksMin[j] + rng.nextInt(tasksMax[i] - tasksMin[i] + 1);
            }
        }
    }

    private StochasticScheduling(
        LocalSolver localsolver,
        int nbTasks,
        int nbMachines,
        int nbScenarios,
        int rngSeed
    ) {
        this.localsolver = localsolver;
        this.nbTasks = nbTasks;
        this.nbMachines = nbMachines;
        this.nbScenarios = nbScenarios;
        generateScenarios(rngSeed);
    }

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

        machineTasks = new LSExpression[nbMachines];
        scenarioMakespan = new LSExpression[nbScenarios];

        // Set decisions: machineTasks[k] represents the tasks processed my machine k
        for (int k = 0; k < nbMachines; k++) {
            machineTasks[k] = model.setVar(nbTasks);
        }

        // Each task must be processed by exactly one machine
        model.constraint(model.partition(machineTasks));

        // Compute makespan of each scenario
        for (int m = 0; m < nbScenarios; m++) {
            LSExpression scenario = model.array(scenarioTaskDurations[m]);
            LSExpression taskSelector = model.lambdaFunction(i -> model.at(scenario, i));
            LSExpression[] machinesMakespan = new LSExpression[nbMachines];

            for (int k = 0; k < nbMachines; k++) {
                machinesMakespan[k] = model.sum(machineTasks[k], taskSelector);
            }
            scenarioMakespan[m] = model.max(machinesMakespan);
        }
        
        // Compute the 90th percentile of scenario makespans.
        LSExpression scenarioMakespanArray = model.array(scenarioMakespan);
        LSExpression sortedScenarioMakespan = model.sort(scenarioMakespanArray);
        stochasticMakespan = model.at(sortedScenarioMakespan, (int)Math.ceil(0.9 * (nbScenarios - 1)));

        model.minimize(stochasticMakespan);
        model.close();

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

        localsolver.solve();
    }

    private void writeSolution() {
        System.out.println();
        System.out.println("Scenario task durations:");
        for (int i = 0; i < nbScenarios; i++) {
            System.out.print("" + i + ": [");
            for (int j = 0; j < nbTasks; ++j) {
                System.out.print("" + scenarioTaskDurations[i][j] + (j == nbTasks - 1 ? "" : ", "));
            }
            System.out.println("]");
        }


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

        for (int m = 0; m < nbMachines; m++) {
            System.out.print("" + m + ": { ");
            LSCollection tasks = machineTasks[m].getCollectionValue();
            for (int i = 0; i < tasks.count(); i++) {
                System.out.print("" + tasks.get(i) + (i == tasks.count() - 1 ? " " : ", "));
            }
            System.out.println("}");
        }
    }

    public static void main(String[] args) {
        try (LocalSolver localsolver = new LocalSolver()) {
            int nbTasks = 10;
            int nbMachines = 2;
            int nbScenarios = 3;
            int rngSeed = 42;
            int timeLimit = 2;
            
            StochasticScheduling model = new StochasticScheduling(
                localsolver,
                nbTasks,
                nbMachines,
                nbScenarios,
                rngSeed
            );

            model.solve(timeLimit);
            model.writeSolution();
        } catch(Exception ex) {
            System.err.println(ex);
            ex.printStackTrace();
            System.exit(1);
        }
    }
};
Compilation/Execution (Windows)
copy %LS_HOME%\bin\*net.dll .
csc StochasticScheduling.cs /reference:localsolvernet.dll
StochasticScheduling
/********** StochasticScheduling.cs **********/

using System;

using localsolver;

public class StochasticScheduling : IDisposable {
    //Number of tasks
    int nbTasks;
   
    //Number of machines
    int nbMachines;
   
    //Number of scenarios
    int nbScenarios;

    //For each scenario, the processing time of each task
    int[][] scenarioTaskDurations;

    // Solver
    LocalSolver localsolver;

    // Decision variable for the assignment of tasks
    LSExpression[] machineTasks;
    
    // For each scenario, the corresponding makespan
    LSExpression[] scenarioMakespan;
    
    // Objective = minimize the 9th decile of all possible makespans
    LSExpression stochasticMakespan;

    private void generateScenarios(int rngSeed) {
        Random rng = new Random(rngSeed);
        
        // Pick random parameters for each task distribution
        int[] tasksMin = new int[nbTasks];
        int[] tasksMax = new int[nbTasks];
        for (int i = 0; i < nbTasks; i++) {
            tasksMin[i] = rng.Next(10, 101);
            tasksMax[i] = tasksMin[i] + rng.Next(51);
        }

        // Sample the distributions to generate the scenarios
        scenarioTaskDurations = new int[nbScenarios][];
        for (int i = 0; i < nbScenarios; i++) {
            scenarioTaskDurations[i] = new int[nbTasks];
            for (int j = 0; j < nbTasks; ++j) {
                scenarioTaskDurations[i][j] = rng.Next(tasksMin[i], tasksMax[i] + 1);
            }
        }
    }

    public StochasticScheduling(
        int nbTasks,
        int nbMachines,
        int nbScenarios,
        int rngSeed
    ) {
        localsolver = new LocalSolver();
        this.nbTasks = nbTasks;
        this.nbMachines = nbMachines;
        this.nbScenarios = nbScenarios;
        generateScenarios(rngSeed);
    }

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

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

        machineTasks = new LSExpression[nbMachines];
        scenarioMakespan = new LSExpression[nbScenarios];

        // Set decisions: machineTasks[k] represents the tasks processed my machine k
        for (int k = 0; k < nbMachines; k++) {
            machineTasks[k] = model.Set(nbTasks);
        }

        // Each task must be processed by exactly one machine
        model.Constraint(model.Partition(machineTasks));

        // Compute makespan of each scenario
        for (int m = 0; m < nbScenarios; m++) {
            LSExpression scenario = model.Array(scenarioTaskDurations[m]);
            LSExpression taskSelector = model.LambdaFunction(i => scenario[i]);
            LSExpression[] machinesMakespan = new LSExpression[nbMachines];

            for (int k = 0; k < nbMachines; k++) {
                machinesMakespan[k] = model.Sum(machineTasks[k], taskSelector);
            }
            scenarioMakespan[m] = model.Max(machinesMakespan);
        }

        // Compute the 90th percentile of scenario makespans.
        LSExpression scenarioMakespanArray = model.Array(scenarioMakespan);
        LSExpression sortedScenarioMakespan = model.Sort(scenarioMakespanArray);
        stochasticMakespan = sortedScenarioMakespan[(int)Math.Ceiling(0.9 * (nbScenarios - 1))];

        model.Minimize(stochasticMakespan);
        model.Close();

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

        localsolver.Solve();
    }

    private void WriteSolution() {
        Console.WriteLine();
        Console.WriteLine("Scenario task durations:");
        for (int i = 0; i < nbScenarios; i++) {
            Console.Write(i + ": [");
            for (int j = 0; j < nbTasks; ++j) {
                Console.Write(scenarioTaskDurations[i][j] + (j == nbTasks - 1 ? "" : ", "));
            }
            Console.WriteLine("]");
        }

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

        for (int m = 0; m < nbMachines; m++) {
            Console.Write(m + ": { ");
            LSCollection tasks = machineTasks[m].GetCollectionValue();
            for (int i = 0; i < tasks.Count(); i++) {
                Console.Write(tasks.Get(i) + (i == tasks.Count() - 1 ? " " : ", "));
            }
            Console.WriteLine("}");
        }
    }

    public static void Main(string[] args) {
        int nbTasks = 10;
        int nbMachines = 2;
        int nbScenarios = 3;
        int rngSeed = 43;
        int timeLimit = 2;
        
        using (StochasticScheduling model = new StochasticScheduling(
            nbTasks,
            nbMachines,
            nbScenarios,
            rngSeed
        )) {
            model.Solve(timeLimit);
            model.WriteSolution();
        }
    }
}