Stochastic Job Shop Scheduling Problem

Problem


In the Job Shop Scheduling Problem, a set of jobs has to be processed on every machine in the shop. The jobs consist in ordered sequences of tasks (called activities). An activity represents the processing of the job on one of the machines and has a given processing time. Each job has one activity per machine, and each activity can only start when the previous one has ended. Each machine can only process one activity at a time. In this example, we consider a stochastic version of the Job Shop Scheduling Problem with multiple scenarios: the processing time of each activity varies from one scenario to another. For a given sequence of jobs and a given scenario, the makespan is the time when all jobs have been processed. The goal of the problem is to find a sequence of jobs that minimizes the maximum makespan over all the scenarios.

Principles learned

Data

The format of the data files is as follows:

  • First line: number of jobs, number of machines, number of scenarios.
  • From the fourth line, for each scenario:
    • For each job: the processing time on each machine (given in processing order).
  • For each job: the processing order (ordered list of visited machines, same for all scenarios).

Program

The Hexaly model for the Stochastic Job Shop Scheduling Problem uses interval decision variables to model the time ranges of the activities in each scenario. We constrain the length of each interval to correspond to the activity’s processing time in the corresponding scenario. We can then write the precedence constraints: for each job and each scenario, each activity of this job must start after the end of the activity processed by the previous machine.

In addition to the interval decisions representing the time ranges of the activities, we also use list decision variables. A list models the ordering of activities on a machine. Using the ‘count’ operator to constrain the size of the lists, we ensure that each machine processes each job. The disjunctive resource constraints — each machine can only process one activity at a time — can be formulated as follows: for all i, the activity processed in position i+1 must start after the end of the activity processed in position i. To model these constraints, we pair up the interval decisions (the time ranges) with the list decisions (the job orderings). We write a lambda function expressing the relationship between two consecutive activities. This function is used within a variadic ‘and’ operator over all activities processed by each machine in each scenario.

The objective is to minimize the maximum makepan among all scenarios, which is the time when all the activities have been processed for each scenario.

Execution
hexaly stochastic_jobshop.hxm inFileName=instances/ft20_10.txt [hxTimeLimit=] [solFileName=]
use io;

/* Read instance data. */
function input() {
    local usage = "Usage: hexaly stochastic_jobshop.hxm inFileName=instanceFile "
            + "[outFileName=outputFile] [hxTimeLimit=timeLimit]";
    if (inFileName == nil) throw usage;

    inFile = io.openRead(inFileName);
    inFile.readln();
    nbJobs = inFile.readInt();
    nbMachines = inFile.readInt();
    nbScenarios = inFile.readInt();
    inFile.readln();
    
    // Processing times for each job on each machine (given in the processing order)
    for [s in 0...nbScenarios][j in 0...nbJobs][m in 0...nbMachines] {
        processingTimesInProcessingOrderPerScenario[s][j][m] = inFile.readInt();
    }

    inFile.readln();
    for [j in 0...nbJobs][k in 0...nbMachines] {
        local m = inFile.readInt()-1;

        // Processing order of machines for each job
        machineOrder[j][k] = m;

        for [s in 0...nbScenarios] {
            // Reorder processing times: processingTime[s][j][m] is the processing time of the
            // task of job j that is processed on machine m in the scenario s
            processingTimePerScenario[s][j][m] = processingTimesInProcessingOrderPerScenario[s][j][k];
        }
    }
    inFile.close();

    // Trivial upper bounds for the start times of the tasks
    maxStart = max[s in 0...nbScenarios](sum[j in 0...nbJobs][m in 0...nbMachines](processingTimePerScenario[s][j][m]));
}


/* Declare the optimization model */
function model() {
    // Interval decisions: time range of each task
    // tasks[s][j][m] is the interval of time of the task of job j which is processed
    // on machine m in the scenario s
    tasks[s in 0...nbScenarios][j in 0...nbJobs][m in 0...nbMachines] <- interval(0, maxStart);

    // Task duration constraints
    for [s in 0...nbScenarios][j in 0...nbJobs][m in 0...nbMachines]
        constraint length(tasks[s][j][m]) == processingTimePerScenario[s][j][m];

    // Precedence constraints between the tasks of a job
    for [s in 0...nbScenarios][j in 0...nbJobs][k in 0...nbMachines-1]
        constraint tasks[s][j][machineOrder[j][k]] < tasks[s][j][machineOrder[j][k + 1]];

    // Sequence of tasks on each machine
    jobsOrder[m in 0...nbMachines] <- list(nbJobs);

    for [m in 0...nbMachines] {
        // Each job has a task scheduled on each machine
        constraint count(jobsOrder[m]) == nbJobs;

        // Disjunctive resource constraints between the tasks on a machine
        for [s in 0...nbScenarios] {
            constraint and(0...nbJobs-1,
                    i => tasks[s][jobsOrder[m][i]][m] < tasks[s][jobsOrder[m][i + 1]][m]);
        }
    }

    // Minimize the maximum makespan: end of the last task of the last job
    // over all scenarios
    makespans[s in 0...nbScenarios] <- max[j in 0...nbJobs](end(tasks[s][j][machineOrder[j][nbMachines - 1]]));
    maxMakespan <- max[s in 0...nbScenarios](makespans[s]);
    minimize maxMakespan;
}

/* Parameterize the solver */
function param() {
    if (hxTimeLimit == nil) hxTimeLimit = 60;
}

/* Write the solution in a file with the following format:
 *  - for each machine, the job sequence */
function output() {
    if (outFileName != nil) {
        outFile = io.openWrite(outFileName);
        println("Solution written in file ", outFileName);
        for [m in 0...nbMachines]
            outFile.println[j in 0...nbJobs](jobsOrder[m].value[j], " ");
    }
}
Execution (Windows)
set PYTHONPATH=%HX_HOME%\bin\python
python stochastic_jobshop.py instances\ft20_10.txt
Execution (Linux)
export PYTHONPATH=/opt/hexaly_13_0/bin/python
python stochastic_jobshop.py instances\ft20_10.txt
import hexaly.optimizer
import sys


def read_instance(filename):
    with open(filename) as f:
        lines = f.readlines()

    first_line = lines[1].split()
    # Number of jobs
    nb_jobs = int(first_line[0])
    # Number of machines
    nb_machines = int(first_line[1])
    # Number of scenarios
    nb_scenarios = int(first_line[2])

    # Processing times for each job on each machine (given in the processing order)
    processing_times_in_processing_order_per_scenario = [[[int(lines[s*(nb_jobs+1)+i].split()[j])
                                                           for j in range(nb_machines)]
                                                          for i in range(3, 3 + nb_jobs)]
                                                         for s in range(nb_scenarios)]

    # Processing order of machines for each job
    machine_order = [[int(lines[i].split()[j]) - 1 for j in range(nb_machines)]
                     for i in range(4 + nb_scenarios*(nb_jobs+1), 4 + nb_scenarios*(nb_jobs+1) + nb_jobs)]

    # Reorder processing times: processing_time[s][j][m] is the processing time of the
    # task of job j that is processed on machine m in the scenario s
    processing_time_per_scenario = [[[processing_times_in_processing_order_per_scenario[s][j][machine_order[j].index(m)]
                                      for m in range(nb_machines)]
                                     for j in range(nb_jobs)]
                                    for s in range(nb_scenarios)]

    # Trivial upper bound for the start times of the tasks
    max_start = max([sum(sum(processing_time_per_scenario[s][j])
                    for j in range(nb_jobs)) for s in range(nb_scenarios)])

    return nb_jobs, nb_machines, nb_scenarios, processing_time_per_scenario, machine_order, max_start


def main(instance_file, output_file, time_limit):
    nb_jobs, nb_machines, nb_scenarios, processing_time_per_scenario, machine_order, max_start = read_instance(
        instance_file)

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

        # Interval decisions: time range of each task
        # tasks[s][j][m] is the interval of time of the task of job j which is processed
        # on machine m in the scenario s
        tasks = [[[model.interval(0, max_start) for m in range(nb_machines)]
                  for j in range(nb_jobs)]
                 for s in range(nb_scenarios)]

        # Task duration constraints
        for s in range(nb_scenarios):
            for j in range(nb_jobs):
                for m in range(0, nb_machines):
                    model.constraint(model.length(tasks[s][j][m]) == processing_time_per_scenario[s][j][m])

        # Create an Hexaly array in order to be able to access it with "at" operators
        task_array = model.array(tasks)

        # Precedence constraints between the tasks of a job
        for s in range(nb_scenarios):
            for j in range(nb_jobs):
                for k in range(nb_machines - 1):
                    model.constraint(
                        tasks[s][j][machine_order[j][k]] < tasks[s][j][machine_order[j][k + 1]])

        # Sequence of tasks on each machine
        jobs_order = [model.list(nb_jobs) for m in range(nb_machines)]

        for m in range(nb_machines):
            # Each job has a task scheduled on each machine
            sequence = jobs_order[m]
            model.constraint(model.eq(model.count(sequence), nb_jobs))

            # Disjunctive resource constraints between the tasks on a machine
            for s in range(nb_scenarios):
                sequence_lambda = model.lambda_function(
                    lambda i: model.lt(model.at(task_array, s, sequence[i], m),
                                       model.at(task_array, s, sequence[i + 1], m)))
                model.constraint(model.and_(model.range(0, nb_jobs - 1), sequence_lambda))

        # Minimize the maximum makespan: end of the last task of the last job
        # over all scenarios
        makespans = [model.max([model.end(tasks[s][j][machine_order[j][nb_machines - 1]]) for j in range(nb_jobs)])
                     for s in range(nb_scenarios)]
        max_makespan = model.max(makespans)
        model.minimize(max_makespan)

        model.close()

        # Parameterize the optimizer
        optimizer.param.time_limit = time_limit

        optimizer.solve()

        #
        # Write the solution in a file with the following format:
        # - for each machine, the job sequence
        #
        if output_file != None:
            final_jobs_order = [list(jobs_order[m].value) for m in range(nb_machines)]
            with open(output_file, "w") as f:
                print("Solution written in file ", output_file)
                for m in range(nb_machines):
                    for j in range(nb_jobs):
                        f.write(str(final_jobs_order[m][j]) + " ")
                    f.write("\n")


if __name__ == '__main__':
    if len(sys.argv) < 2:
        print("Usage: python stochastic_jobshop.py instance_file [output_file] [time_limit]")
        sys.exit(1)

    instance_file = sys.argv[1]
    output_file = sys.argv[2] if len(sys.argv) >= 3 else None
    time_limit = int(sys.argv[3]) if len(sys.argv) >= 4 else 60
    main(instance_file, output_file, time_limit)
Compilation / Execution (Windows)
cl /EHsc stochastic_jobshop.cpp -I%HX_HOME%\include /link %HX_HOME%\bin\hexaly130.lib
stochastic_jobshop instances\ft20_10.txt
Compilation / Execution (Linux)
g++ stochastic_jobshop.cpp -I/opt/hexaly_13_0/include -lhexaly130 -lpthread -o stochastic_jobshop
./stochastic_jobshop instances\ft20_10.txt
#include "optimizer/hexalyoptimizer.h"
#include <algorithm>
#include <fstream>
#include <iostream>
#include <limits>
#include <numeric>
#include <vector>

using namespace hexaly;
using namespace std;

class StochasticJobshop {
private:
    // Number of jobs
    int nbJobs;
    // Number of machines
    int nbMachines;
    // Number of scenarios
    int nbScenarios;
    // Processing order of machines for each job
    vector<vector<int>> machineOrder;
    // Processing time on each machine for each job (given in the machine order and for a given scenario)
    vector<vector<vector<int>>> processingTimePerScenario;
    // Trivial upper bound for the start times of the tasks
    int maxStart;

    // Hexaly Optimizer
    HexalyOptimizer optimizer;
    // Decision variables: time range of each task
    vector<vector<vector<HxExpression>>> tasks;
    // Decision variables: sequence of tasks on each machine
    vector<HxExpression> jobsOrder;
    // Objective = minimize the maximum of all makespans
    HxExpression maxMakespan;

public:
    StochasticJobshop() : optimizer() {}

    void readInstance(const string& fileName) {

        ifstream infile;
        infile.exceptions(ifstream::failbit | ifstream::badbit);
        infile.open(fileName.c_str());

        infile.ignore(numeric_limits<streamsize>::max(), '\n');
        infile >> nbJobs;
        infile >> nbMachines;
        infile >> nbScenarios;
        infile.ignore(numeric_limits<streamsize>::max(), '\n');

        // Processing times for each job on each machine (given in the processing order)
        infile.ignore(numeric_limits<streamsize>::max(), '\n');
        vector<vector<vector<int>>> processingTimeInProcessingOrderPerScenario =
            vector<vector<vector<int>>>(nbScenarios, vector<vector<int>>(nbJobs, vector<int>(nbMachines)));
            
        for (int s = 0; s < nbScenarios; ++s) {
            for (int j = 0; j < nbJobs; ++j) {
                for (int m = 0; m < nbMachines; ++m) {
                    infile >> processingTimeInProcessingOrderPerScenario[s][j][m];
                }
            }
            infile.ignore(numeric_limits<streamsize>::max(), '\n');
        }

        // Processing order of machines for each job
        infile.ignore(numeric_limits<streamsize>::max(), '\n');
        infile.ignore(numeric_limits<streamsize>::max(), '\n');
        machineOrder.resize(nbJobs);
        for (int j = 0; j < nbJobs; ++j) {
            machineOrder[j].resize(nbMachines);
            for (int m = 0; m < nbMachines; ++m) {
                int x;
                infile >> x;
                machineOrder[j][m] = x - 1;
            }
        }

        // Reorder processing times: processingTimePerScenario[s][j][m] is the processing time of the
        // task of job j that is processed on machine m in scenario s
        for (int s = 0; s < nbScenarios; ++s) {
            processingTimePerScenario.resize(nbScenarios);
            for (int j = 0; j < nbJobs; ++j) {
                processingTimePerScenario[s].resize(nbJobs);
                for (int m = 0; m < nbMachines; ++m) {
                    processingTimePerScenario[s][j].resize(nbMachines);
                    vector<int>::iterator findM = find(machineOrder[j].begin(), machineOrder[j].end(), m);
                    unsigned int k = distance(machineOrder[j].begin(), findM);
                    processingTimePerScenario[s][j][m] = processingTimeInProcessingOrderPerScenario[s][j][k];
                }
            }
        }

        // Trivial upper bound for the start times of the tasks
        vector<int> maxStartPerScenario = vector<int>(nbScenarios);
        for (int s = 0; s < nbScenarios; ++s) {
            for (int j = 0; j < nbJobs; ++j) {
                maxStartPerScenario[s] +=
                    accumulate(processingTimePerScenario[s][j].begin(), processingTimePerScenario[s][j].end(), 0);
            }
        }
        maxStart = *max_element(maxStartPerScenario.begin(), maxStartPerScenario.end());
        infile.close();
    }
    
    void solve(int timeLimit) {
        // Declare the optimization model
        HxModel model = optimizer.getModel();

        // Interval decisions: time range of each task
        // tasks[s][j][m] is the interval of time of the task of job j which is 
        // processed on machine m in scenario s
        tasks.resize(nbScenarios);
        for (unsigned int s = 0; s < nbScenarios; ++s) {
            tasks[s].resize(nbJobs);
            for (unsigned int j = 0; j < nbJobs; ++j) {
                tasks[s][j].resize(nbMachines);
                for (unsigned int m = 0; m < nbMachines; ++m) {
                    tasks[s][j][m] = model.intervalVar(0, maxStart);

                    // Task duration constraints
                    model.constraint(model.length(tasks[s][j][m]) == processingTimePerScenario[s][j][m]);
                }
            }
        }

        // Create an Hexaly array in order to be able to access it with "at" operators
        vector<HxExpression> taskArray = vector<HxExpression>(nbScenarios);
        for (int s = 0; s < nbScenarios; ++s) {
            taskArray[s] = model.array();
            for (int j = 0; j < nbJobs; ++j) {
                taskArray[s].addOperand(model.array(tasks[s][j].begin(), tasks[s][j].end()));
            }
        }

        // Precedence constraints between the tasks of a job
        for (int s = 0; s < nbScenarios; ++s) {
            for (int j = 0; j < nbJobs; ++j) {
                for (int k = 0; k < nbMachines - 1; ++k) {
                    model.constraint(tasks[s][j][machineOrder[j][k]] < tasks[s][j][machineOrder[j][k + 1]]);
                }
            }
        }

        // Sequence of tasks on each machine
        jobsOrder.resize(nbMachines);
        for (int m = 0; m < nbMachines; ++m) {
            jobsOrder[m] = model.listVar(nbJobs);
        }

        for (int m = 0; m < nbMachines; ++m) {
            // Each job has a task scheduled on each machine
            HxExpression sequence = jobsOrder[m];
            model.constraint(model.eq(model.count(sequence), nbJobs));

            // Disjunctive resource constraints between the tasks on a machine
            for (int s = 0; s < nbScenarios; ++s) {
                HxExpression sequenceLambda = model.createLambdaFunction([&](HxExpression i) {
                    return model.at(taskArray[s], sequence[i], m) < model.at(taskArray[s], sequence[i + 1], m);
                });
                model.constraint(model.and_(model.range(0, nbJobs - 1), sequenceLambda));
            }
        }

        // Minimize the maximum makespan: end of the last task of the last job
        // over all scenarios
        vector<HxExpression> makespans = vector<HxExpression>(nbScenarios);
        for (int s = 0; s < nbScenarios; ++s) {
            makespans[s] = model.max();
            for (int j = 0; j < nbJobs; ++j) {
                makespans[s].addOperand(model.end(tasks[s][j][machineOrder[j][nbMachines - 1]]));
            }
        }
        maxMakespan = model.max();
        for (int s = 0; s < nbScenarios; ++s) 
            maxMakespan.addOperand(makespans[s]);

        model.minimize(maxMakespan);
        model.close();

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

        optimizer.solve();
    }

    /* Write the solution in a file with the following format:
     *  - for each machine, the job sequence */
    void writeSolution(const string& fileName) {
        ofstream outfile;
        outfile.exceptions(ofstream::failbit | ofstream::badbit);
        outfile.open(fileName.c_str());
        cout << "Solution written in file " << fileName << endl;

        for (int m = 0; m < nbMachines; ++m) {
            HxCollection finalJobsOrder = jobsOrder[m].getCollectionValue();
            for (int j = 0; j < nbJobs; ++j) {
                outfile << finalJobsOrder.get(j) << " ";
            }
            outfile << endl;
        }
        outfile.close();
    }
};

int main(int argc, char** argv) {
    if (argc < 2) {
        cout << "Usage: stochastic_jobshop instanceFile [outputFile] [timeLimit]" << endl;
        exit(1);
    }

    const char* instanceFile = argv[1];
    const char* outputFile = argc > 2 ? argv[2] : NULL;
    const char* strTimeLimit = argc > 3 ? argv[3] : "60";

    StochasticJobshop model;
    try {
        model.readInstance(instanceFile);
        const int timeLimit = atoi(strTimeLimit);
        model.solve(timeLimit);
        if (outputFile != NULL)
            model.writeSolution(outputFile);
        return 0;
    } catch (const exception& e) {
        cerr << "An error occurred: " << e.what() << endl;
        return 1;
    }
}
Compilation / Execution (Windows)
copy %HX_HOME%\bin\Hexaly.NET.dll .
csc StochasticJobshop.cs /reference:Hexaly.NET.dll
StochasticJobshop instances\ft20_10.txt
using System;
using System.Linq;
using System.IO;
using Hexaly.Optimizer;

public class StochasticJobshop : IDisposable
{
    // Number of jobs
    private int nbJobs;

    // Number of machines
    private int nbMachines;

    // Number of machines
    private int nbScenarios;

    // Processing order of machines for each job
    private int[,] machineOrder;

    // Processing time on each machine for each job (given in the machine order)
    private long[,,] processingTimePerScenario;

    // Trivial upper bound for the start times of the tasks
    private long maxStart;

    // Hexaly Optimizer
    private HexalyOptimizer optimizer;

    // Decision variables: time range of each task
    private HxExpression[,,] tasks;

    // Decision variables: sequence of tasks on each machine
    private HxExpression[] jobsOrder;

    // Objective = minimize the maximum of all makespans
    private HxExpression maxMakespan;

    public StochasticJobshop()
    {
        optimizer = new HexalyOptimizer();
    }

    public void ReadInstance(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            input.ReadLine();
            string[] splitted = input.ReadLine().Split(' ');
            nbJobs = int.Parse(splitted[0]);
            nbMachines = int.Parse(splitted[1]);
            nbScenarios = int.Parse(splitted[2]);

            // Processing times for each job on each machine (given in the processing order)
            input.ReadLine();
            long[,,] processingTimeInProcessingOrderPerScenario = new long[
                nbScenarios,
                nbJobs,
                nbMachines
            ];
            for (int s = 0; s < nbScenarios; ++s)
            {
                for (int j = 0; j < nbJobs; ++j)
                {
                    splitted = input.ReadLine().Trim().Split(' ');
                    for (int m = 0; m < nbMachines; ++m)
                        processingTimeInProcessingOrderPerScenario[s, j, m] = long.Parse(
                            splitted[m]
                        );
                }
                input.ReadLine();
            }

            // Processing order of machines for each job
            input.ReadLine();
            machineOrder = new int[nbJobs, nbMachines];
            for (int j = 0; j < nbJobs; ++j)
            {
                splitted = input.ReadLine().Trim().Split(' ');
                for (int m = 0; m < nbMachines; ++m)
                    machineOrder[j, m] = int.Parse(splitted[m]) - 1;
            }

            // Reorder processing times: processingTimePerScenario[s, j, m] is the processing time of the
            // task of job j that is processed on machine m in scenario s
            processingTimePerScenario = new long[nbScenarios, nbJobs, nbMachines];
            // Trivial upper bound for the start times of the tasks
            long[] maxStartPerScenario = new long[nbScenarios];
            for (int s = 0; s < nbScenarios; ++s)
            {
                for (int j = 0; j < nbJobs; ++j)
                {
                    for (int m = 0; m < nbMachines; ++m)
                    {
                        int machineIndex = nbMachines;
                        for (int k = 0; k < nbMachines; ++k)
                        {
                            if (machineOrder[j, k] == m)
                            {
                                machineIndex = k;
                                break;
                            }
                        }
                        processingTimePerScenario[s, j, m] =
                            processingTimeInProcessingOrderPerScenario[s, j, machineIndex];
                        maxStartPerScenario[s] += processingTimePerScenario[s, j, m];
                    }
                }
            }
            maxStart = maxStartPerScenario.Max();
        }
    }

    public void Dispose()
    {
        optimizer.Dispose();
    }

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

        // Interval decisions: time range of each task
        // tasks[s][j][m] is the interval of time of the task of job j which is processed on machine m
        // in scenario s
        tasks = new HxExpression[nbScenarios, nbJobs, nbMachines];

        for (int s = 0; s < nbScenarios; ++s)
        {
            for (int j = 0; j < nbJobs; ++j)
            {
                for (int m = 0; m < nbMachines; ++m)
                {
                    tasks[s, j, m] = model.Interval(0, maxStart);

                    // Task duration constraints
                    model.Constraint(
                        model.Length(tasks[s, j, m]) == processingTimePerScenario[s, j, m]
                    );
                }
            }
        }

        // Create a HexalyOptimizer array in order to be able to access it with "at" operators
        HxExpression taskArray = model.Array(tasks);

        // Precedence constraints between the tasks of a job
        for (int s = 0; s < nbScenarios; ++s)
        {
            for (int j = 0; j < nbJobs; ++j)
            {
                for (int k = 0; k < nbMachines - 1; ++k)
                {
                    model.Constraint(
                        tasks[s, j, machineOrder[j, k]] < tasks[s, j, machineOrder[j, k + 1]]
                    );
                }
            }
        }

        // Sequence of tasks on each machine
        jobsOrder = new HxExpression[nbMachines];
        for (int m = 0; m < nbMachines; ++m)
            jobsOrder[m] = model.List(nbJobs);

        for (int m = 0; m < nbMachines; ++m)
        {
            // Each job has a task scheduled on each machine
            HxExpression sequence = jobsOrder[m];
            model.Constraint(model.Count(sequence) == nbJobs);

            // Disjunctive resource constraints between the tasks on a machine
            for (int s = 0; s < nbScenarios; ++s)
            {
                HxExpression sequenceLambda = model.LambdaFunction(
                    i => taskArray[s][sequence[i], m] < taskArray[s][sequence[i + 1], m]
                );
                model.Constraint(model.And(model.Range(0, nbJobs - 1), sequenceLambda));
            }
        }

        // Minimize the maximum makespan: end of the last task of the last job
        // over all scenarios
        HxExpression[] makespans = new HxExpression[nbScenarios];
        for (int s = 0; s < nbScenarios; ++s)
        {
            makespans[s] = model.Max();
            for (int j = 0; j < nbJobs; ++j)
            {
                makespans[s].AddOperand(model.End(tasks[s, j, machineOrder[j, nbMachines - 1]]));
            }
        }

        maxMakespan = model.Max();
        for (int s = 0; s < nbScenarios; ++s)
            maxMakespan.AddOperand(makespans[s]);

        model.Minimize(maxMakespan);
        model.Close();

        // Parameterize the optimizer
        optimizer.GetParam().SetTimeLimit(timeLimit);

        optimizer.Solve();
    }

    /* Write the solution in a file with the following format:
     *  - for each machine, the job sequence */
    public void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            Console.WriteLine("Solution written in file " + fileName);
            for (int m = 0; m < nbMachines; ++m)
            {
                HxCollection finalJobsOrder = jobsOrder[m].GetCollectionValue();
                for (int i = 0; i < nbJobs; ++i)
                {
                    int j = (int)finalJobsOrder.Get(i);
                    output.Write(j + " ");
                }
                output.WriteLine();
            }
        }
    }

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: StochasticJobshop instanceFile [outputFile] [timeLimit]");
            System.Environment.Exit(1);
        }

        string instanceFile = args[0];
        string outputFile = args.Length > 1 ? args[1] : null;
        string strTimeLimit = args.Length > 2 ? args[2] : "60";

        using (StochasticJobshop model = new StochasticJobshop())
        {
            model.ReadInstance(instanceFile);
            model.Solve(int.Parse(strTimeLimit));
            if (outputFile != null)
                model.WriteSolution(outputFile);
        }
    }
}
Compilation / Execution (Windows)
javac StochasticJobshop.java -cp %HX_HOME%\bin\hexaly.jar
java -cp %HX_HOME%\bin\hexaly.jar;. StochasticJobshop instances\ft20_10.txt
Compilation / Execution (Linux)
javac StochasticJobshop.java -cp /opt/hexaly_13_0/bin/hexaly.jar
java -cp /opt/hexaly_13_0/bin/hexaly.jar:. StochasticJobshop instances\ft20_10.txt
import java.util.*;
import java.io.*;
import com.hexaly.optimizer.*;

public class StochasticJobshop {
    // Number of jobs
    private int nbJobs;
    // Number of machines
    private int nbMachines;
    // Number of scenarios
    private int nbScenarios;
    // Processing time on each machine for each job (given in the machine order)
    private long[][][] processingTimePerScenario;
    // Processing order of machines for each job
    private int[][] machineOrder;
    // Trivial upper bound for the start times of the tasks
    private long maxStart;

    // Hexaly Optimizer
    final HexalyOptimizer optimizer;
    // Decision variables: time range of each task
    private HxExpression[][][] tasks;
    // Decision variables: sequence of tasks on each machine
    private HxExpression[] jobsOrder;
    // Objective = minimize the maximum of all makespans
    private HxExpression maxMakespan;

    public StochasticJobshop(HexalyOptimizer optimizer) throws IOException {
        this.optimizer = optimizer;
    }

    public void readInstance(String fileName) throws IOException {
        try (Scanner input = new Scanner(new File(fileName))) {
            input.useLocale(Locale.ROOT);
            input.nextLine();
            nbJobs = input.nextInt();
            nbMachines = input.nextInt();
            nbScenarios = input.nextInt();

            input.nextLine();
            input.nextLine();
            // Processing times for each job on each machine (given in the processing order)
            long[][][] processingTimesInProcessingOrderPerScenario = new long[nbScenarios][nbJobs][nbMachines];
            for (int s = 0; s < nbScenarios; ++s) {
                for (int j = 0; j < nbJobs; ++j) {
                    for (int m = 0; m < nbMachines; ++m) {
                        processingTimesInProcessingOrderPerScenario[s][j][m] = input.nextInt();
                    }
                }
                input.nextLine();
            }
            // Processing order of machines for each job
            input.nextLine();
            input.nextLine();
            machineOrder = new int[nbJobs][nbMachines];
            for (int j = 0; j < nbJobs; ++j) {
                for (int m = 0; m < nbMachines; ++m) {
                    machineOrder[j][m] = input.nextInt() - 1;
                }
            }

            // Reorder processing times: processingTimePerScenario[s][j][m] is the
            // processing time of the task of job j that is processed on machine m for a
            // given scenario s
            processingTimePerScenario = new long[nbScenarios][nbJobs][nbMachines];

            // Trivial upper bound for the start times of the tasks
            long[] maxStartPerScenario = new long[nbScenarios];
            for (int s = 0; s < nbScenarios; ++s) {
                maxStartPerScenario[s] = 0;
                for (int j = 0; j < nbJobs; ++j) {
                    for (int m = 0; m < nbMachines; ++m) {
                        int machineIndex = nbMachines;
                        for (int k = 0; k < nbMachines; ++k) {
                            if (machineOrder[j][k] == m) {
                                machineIndex = k;
                                break;
                            }
                        }
                        processingTimePerScenario[s][j][m] = processingTimesInProcessingOrderPerScenario[s][j][machineIndex];
                        maxStartPerScenario[s] += processingTimePerScenario[s][j][m];
                    }
                }
            }
            maxStart = Arrays.stream(maxStartPerScenario).max().getAsLong();
        }
    }

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

        // Interval decisions: time range of each task tasks[s][j][m] is the interval of
        // time of the task of job j which is processed on machine m in the scenario s
        tasks = new HxExpression[nbScenarios][nbJobs][nbMachines];
        for (int j = 0; j < nbJobs; ++j) {
            for (int m = 0; m < nbMachines; ++m) {
                for (int s = 0; s < nbScenarios; ++s) {
                    tasks[s][j][m] = model.intervalVar(0, maxStart);
                    // Task duration constraints
                    model.constraint(model.eq(model.length(tasks[s][j][m]), processingTimePerScenario[s][j][m]));
                }
            }
        }

        // Create a HexalyOptimizer array in order to be able to access it with "at"
        // operators
        HxExpression taskArray = model.array(tasks);

        // Precedence constraints between the tasks of a job
        for (int s = 0; s < nbScenarios; ++s) {
            for (int j = 0; j < nbJobs; ++j) {
                for (int k = 0; k < nbMachines - 1; ++k) {
                    model.constraint(model.lt(tasks[s][j][machineOrder[j][k]], tasks[s][j][machineOrder[j][k + 1]]));
                }
            }
        }

        // Sequence of tasks on each machine
        jobsOrder = new HxExpression[nbMachines];
        for (int m = 0; m < nbMachines; ++m) {
            jobsOrder[m] = model.listVar(nbJobs);
        }

        for (int m = 0; m < nbMachines; ++m) {
            // Each job has a task scheduled on each machine
            HxExpression sequence = jobsOrder[m];
            model.constraint(model.eq(model.count(sequence), nbJobs));

            // Disjunctive resource constraints between the tasks on a machine
            for (int s = 0; s < nbScenarios; ++s) {
                HxExpression mExpr = model.createConstant(m);
                HxExpression sExpr = model.createConstant(s);
                HxExpression sequenceLambda = model
                        .lambdaFunction(i -> model.lt(model.at(taskArray, sExpr, model.at(sequence, i), mExpr),
                                model.at(taskArray, sExpr, model.at(sequence, model.sum(i, 1)), mExpr)));
                model.constraint(model.and(model.range(0, nbJobs - 1), sequenceLambda));
            }
        }

        // Minimize the maximum makespan: end of the last task of the last job
        // over all scenarios
        HxExpression[] makespans = new HxExpression[nbScenarios];
        for (int s = 0; s < nbScenarios; ++s) {
            makespans[s] = model.max();
            for (int j = 0; j < nbJobs; ++j) {
                makespans[s].addOperand(model.end(tasks[s][j][machineOrder[j][nbMachines - 1]]));
            }
        }

        maxMakespan = model.max();
        for (int s = 0; s < nbScenarios; ++s) {
            maxMakespan.addOperand(makespans[s]);
        }

        model.minimize(maxMakespan);

        model.close();

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

        optimizer.solve();
    }

    /*
     * Write the solution in a file with the following format:
     * - for each machine, the job sequence
     */
    public void writeSolution(String fileName) throws IOException {
        try (PrintWriter output = new PrintWriter(fileName)) {
            System.out.println("Solution written in file " + fileName);

            for (int m = 0; m < nbMachines; ++m) {
                HxCollection finalJobsOrder = jobsOrder[m].getCollectionValue();
                for (int i = 0; i < nbJobs; ++i) {
                    int j = Math.toIntExact(finalJobsOrder.get(i));
                    output.write(j + " ");
                }
                output.write("\n");
            }
        }
    }

    public static void main(String[] args) {
        if (args.length < 1) {
            System.out.println("Usage: java StochasticJobshop instanceFile [outputFile] [timeLimit]");
            System.exit(1);
        }

        String instanceFile = args[0];
        String outputFile = args.length > 1 ? args[1] : null;
        String strTimeLimit = args.length > 2 ? args[2] : "60";

        try (HexalyOptimizer optimizer = new HexalyOptimizer()) {
            StochasticJobshop model = new StochasticJobshop(optimizer);
            model.readInstance(instanceFile);
            model.solve(Integer.parseInt(strTimeLimit));
            if (outputFile != null) {
                model.writeSolution(outputFile);
            }
        } catch (Exception ex) {
            System.err.println(ex);
            ex.printStackTrace();
            System.exit(1);
        }
    }
}