Open Shop Scheduling Problem

Principles learned

  • Use lambda functions in standard methods

  • Use variadic operators

  • Use interval decision variables

  • Order interval decision variables by pairing them up with a list variable

Problem

../_images/openshop-scheduling-problem.svg

In the open shop scheduling problem, a set of jobs has to be processed on every machine of the shop. Each job consists in an unordered sequence of tasks (called activities), each representing the processing of the job on one of the machines. Each job has one activity per machine, and cannot start an activity while another activity of the job is running. Each activity has a given processing time and each machine can only process one activity at a time.

The goal is to find a sequence of jobs that minimizes the makespan: the time when all jobs have been processed.

Download the example


Data

The instances provided are from Taillard. The format of the data files is as follows:

  • First line: number of jobs, number of machines, seed used to generate the instance, upper and lower bound previously found.

  • For each job: the processing time of each activity on its assigned machine

  • For each job: the machine id assigned to each activity

Program

We use interval decision variables to model the time ranges of the activities. The length of the interval is constrained by the processing time of each activity.

In addition to the interval decisions representing the time ranges of the activities, we also use list decision variables. As in the Jobshop example, a list models the ordering of activities within a machine. We constrain all the jobs to be processed on each machine thanks to the count operator.

The disjunctive resource constraints — each machine can only process one activity at a time — can be reformulated as follows: given a sequence of jobs, the activity corresponding to any job must be placed after the activity corresponding to the previous job.

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 the two consecutive activities. This function is used within an and operator over all activities processed by a machine.

To model the disjunctive activity constraints - each job can only have one running activity at a time - we use the same strategy: a list models the ordering of activities within a job. We constrain all the activities to be processed thanks to the “count” operator. Similarly to the disjunctive resource constraints, we use an “and” operator over all machines for a given job. Hence, for a given job only one activity is executed at a time.

The makespan to minimize is the time when all the activities have been processed; the maximum end time for all machines.

If you are interested in more specific cases, where the activities are ordered on each machine, you can now study our jobshop model or our flowshop model.

Execution:
hexaly openshop.hxm inFileName=instances/tai2020_5.txt [hxTimeLimit=] [solFileName=]
/********** Openshop **********/

use io;

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

    inFile = io.openRead(inFileName);
    inFile.readln();
    nbJobs = inFile.readInt();
    nbMachines = inFile.readInt();
    inFile.readln();
    inFile.readln();
    // Processing times for each job on each machine
    // (given in the task order, the processing order is a decision variable)
    processingTimesActivityOrder[j in 0...nbJobs][m in 0...nbMachines] = inFile.readInt();
    inFile.readln();
    // Reorder processing times: processingTime[j][m] is the processing time of the
    // task of job j that is processed on machine m
    for [j in 0...nbJobs][k in 0...nbMachines] {
        local m = inFile.readInt() - 1;
        processingTimes[j][m] = processingTimesActivityOrder[j][k];
    }

    inFile.close();

    maxStart = sum[j in 0...nbJobs][m in 0...nbMachines](processingTimes[j][m]);
}

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

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

    // List of the jobs on each machine
    jobsOrder[m in 0...nbMachines] <- list(nbJobs);
    for [m in 0...nbMachines] {
        // Each job is scheduled on every machine
        constraint count(jobsOrder[m]) == nbJobs;

        // Every machine executes a single task at a time
        constraint and(0...nbJobs - 1, i => tasks[jobsOrder[m][i]][m] < tasks[jobsOrder[m][i + 1]][m]);
    }

    // List of the machines for each job task
    machinesOrder[j in 0...nbJobs] <- list(nbMachines);
    for [j in 0...nbJobs] {
        // Every task is scheduled on its corresponding machine
        constraint count(machinesOrder[j]) == nbMachines;

        // A job has a single task at a time
        constraint and(0...nbMachines - 1, k => tasks[j][machinesOrder[j][k]] < tasks[j][machinesOrder[j][k + 1]]);
    }

    // Minimize the makespan: the end of the last task
    makespan <- max[m in 0...nbMachines][j in 0...nbJobs](end(tasks[j][m]));
    minimize makespan;
}

// Parametrize the solver
function param() {
    if (hxTimeLimit == nil) hxTimeLimit = 5;
}

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


def read_instance(filename):
    # The input files follow the "Taillard" format
    with open(filename, 'r') as f:
        lines = f.readlines()

    first_line = lines[1].split()
    nb_jobs = int(first_line[0])
    nb_machines = int(first_line[1])

    # Processing times for each job on each machine
    # (given in the task order, the processing order is a decision variable)
    processing_times_task_order = [[int(proc_time) for proc_time in line.split()]
                                   for line in lines[3:3 + nb_jobs]]

    # Index of machines for each task
    machine_index = [[int(machine_i) - 1 for machine_i in line.split()]
                     for line in lines[4 + nb_jobs:4 + 2 * nb_jobs]]

    # Reorder processing times: processingTime[j][m] is the processing time of the
    # task of job j that is processed on machine m
    processing_times = [[processing_times_task_order[j][machine_index[j].index(m)]
                         for m in range(nb_machines)] for j in range(nb_jobs)]

    # Trivial upper bound for the start time of tasks
    max_start = sum(map(lambda processing_times_job: sum(processing_times_job), processing_times))

    return nb_jobs, nb_machines, processing_times, max_start


def main(instance_file, output_file, time_limit):
    nb_jobs, nb_machines, processing_times, 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[j][m] is the interval of time of the task of job j
        # which is processed on machine m
        tasks = [[model.interval(0, max_start) for _ in range(nb_machines)] for _ in range(nb_jobs)]

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

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

        # List of the jobs on each machine
        jobs_order = [model.list(nb_jobs) for _ in range(nb_machines)]
        for m in range(nb_machines):
            # Each job is scheduled on every machine
            model.constraint(model.eq(model.count(jobs_order[m]), nb_jobs))

            # Every machine executes a single task at a time
            sequence_lambda = model.lambda_function(lambda i:
                model.at(task_array, jobs_order[m][i], m) < model.at(task_array, jobs_order[m][i + 1], m))
            model.constraint(model.and_(model.range(0, nb_jobs - 1), sequence_lambda))

        # List of the machines for each job
        machines_order = [model.list(nb_machines) for _ in range(nb_jobs)]
        for j in range(nb_jobs):
            # Every task is scheduled on its corresponding machine
            model.constraint(model.eq(model.count(machines_order[j]), nb_machines))

            # A job has a single task at a time
            sequence_lambda = model.lambda_function(lambda k:
                    model.at(task_array, j, machines_order[j][k]) < model.at(task_array, j, machines_order[j][k + 1]))
            model.constraint(model.and_(model.range(0, nb_machines - 1), sequence_lambda))

        # Minimize the makespan: the end of the last task
        makespan = model.max([model.end(model.at(task_array, j, m))
                             for j in range(nb_jobs) for m in range(nb_machines)])
        model.minimize(makespan)

        model.close()

        # Parametrize 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 is not None:
            with open(output_file, 'w') as f:
                for m in range(nb_machines):
                    line = ""
                    for j in range(nb_jobs):
                        line += str(jobs_order[m].value[j]) + " "
                    f.write(line + "\n")
            print("Solution written in file ", output_file)


if __name__ == '__main__':
    if len(sys.argv) < 2:
        print(
            "Usage: python openshop.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 openshop.cpp -I%HX_HOME%\include /link %HX_HOME%\bin\hexaly130.lib
openshop instances\tai2020_5.txt
Compilation / Execution (Linux)
g++ openshop.cpp -I/opt/hexaly_13_0/include -lhexaly130 -lpthread -o openshop
./openshop instances/tai2020_5.txt
#include "optimizer/hexalyoptimizer.h"
#include <algorithm>
#include <fstream>
#include <iostream>
#include <limits>
#include <numeric>
#include <vector>

using namespace hexaly;

class Openshop {
private:
    // Number of jobs
    int nbJobs;
    // Number of machines
    int nbMachines;
    // Processing time on each machine for each job task
    std::vector<std::vector<int>> processingTime;
    // Trivial upper bound for the start time of tasks
    int maxStart;

    // Hexaly Optimizer
    HexalyOptimizer optimizer;
    // Decision variables : time range of each task
    std::vector<std::vector<HxExpression>> tasks;
    // Decision variables : processing order of jobs for each machine
    std::vector<HxExpression> jobsOrder;
    // Decision variables : processing order of machines for each job
    std::vector<HxExpression> machinesOrder;
    // Objective = minimize the makespan: end of the last task of the last job
    HxExpression makespan;

public:
    Openshop() : optimizer() {}

    // The input files follow the "Taillard" format
    void readInstance(const std::string& instanceFile) {
        std::ifstream infile;
        infile.exceptions(std::ifstream::failbit | std::ifstream::badbit);
        infile.open(instanceFile.c_str());

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

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

        // Index of machines for each task
        infile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
        std::vector<std::vector<int>> machineIndex =
            std::vector<std::vector<int>>(nbJobs, std::vector<int>(nbMachines));
        for (int j = 0; j < nbJobs; ++j) {
            for (int m = 0; m < nbMachines; ++m) {
                int x;
                infile >> x;
                machineIndex[j][m] = x - 1;
            }
        }

        infile.close();

        // Reorder processing times: processingTime[j][m] is the processing time of the
        // task of job j that is processed on machine m
        processingTime.resize(nbJobs);
        for (int j = 0; j < nbJobs; ++j) {
            processingTime[j].resize(nbMachines);
            for (int m = 0; m < nbMachines; ++m) {
                std::vector<int>::iterator findM = std::find(machineIndex[j].begin(), machineIndex[j].end(), m);
                int k = std::distance(machineIndex[j].begin(), findM);
                processingTime[j][m] = processingTimesActivityOrder[j][k];
            }
        }

        // Trivial upper bound for the start time of tasks
        maxStart = 0;
        for (int j = 0; j < nbJobs; ++j) {
            maxStart += std::accumulate(processingTime[j].begin(), processingTime[j].end(), 0);
        }
    }

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

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

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

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

        // Sequence of tasks on each machine
        jobsOrder.resize(nbMachines);
        for (int m = 0; m < nbMachines; ++m) {
            jobsOrder[m] = model.listVar(nbJobs);
            // Each job is scheduled on every machine
            model.constraint(model.eq(model.count(jobsOrder[m]), nbJobs));

            // Every machine executes a single task at a time
            HxExpression sequenceLambda = model.createLambdaFunction([&](HxExpression i) {
                return model.at(taskArray, jobsOrder[m][i], m) < model.at(taskArray, jobsOrder[m][i + 1], m);
            });
            model.constraint(model.and_(model.range(0, nbJobs - 1), sequenceLambda));
        }

        // Sequence of machines for each job
        machinesOrder.resize(nbJobs);
        for (int j = 0; j < nbJobs; ++j) {
            machinesOrder[j] = model.listVar(nbMachines);
            // Every task is scheduled on its corresponding machine
            model.constraint(model.eq(model.count(machinesOrder[j]), nbMachines));

            // A job has a single task at a time
            HxExpression sequenceLambda = model.createLambdaFunction([&](HxExpression k) {
                return model.at(taskArray, j, machinesOrder[j][k]) < model.at(taskArray, j, machinesOrder[j][k + 1]);
            });
            model.constraint(model.and_(model.range(0, nbMachines - 1), sequenceLambda));
        }

        // Minimize the makespan: the end of the last task
        makespan = model.max();
        for (int m = 0; m < nbMachines; ++m) {
            for (int j = 0; j < nbJobs; ++j) {
                makespan.addOperand(model.end(model.at(taskArray, j, m)));
            }
        }
        model.minimize(makespan);

        model.close();

        // Parametrize 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 std::string& fileName) {
        std::ofstream outfile(fileName.c_str());
        if (!outfile.is_open()) {
            std::cerr << "File " << fileName << " cannot be opened." << std::endl;
            exit(1);
        }
        for (int m = 0; m < nbMachines; ++m) {
            HxCollection finalJobsOrder = jobsOrder[m].getCollectionValue();
            for (int j = 0; j < nbJobs; ++j) {
                outfile << finalJobsOrder.get(j) << " ";
            }
            outfile << std::endl;
        }
        outfile.close();
        std::cout << "Solution written in file " << fileName << std::endl;
    }
};

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

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

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

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

    // Number of machines
    private int nbMachines;

    // Processing time on each machine for each job task
    private long[,] processingTime;

    // 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 : processing order of jobs for each machine
    private HxExpression[] jobsOrder;

    // Decision variables : processing order of machines for each jobs
    private HxExpression[] machinesOrder;

    // Objective = minimize the makespan: end of the last task of the last job
    private HxExpression makespan;

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

    // The input files follow the "Taillard" format
    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]);

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

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

            // Reorder processing times: processingTime[j, m] is the processing time of the
            // task of job j that is processed on machine m
            processingTime = new long[nbJobs, nbMachines];
            // Trivial upper bound for the start times of the tasks
            maxStart = 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 (machineIndexes[j, k] == m)
                        {
                            machineIndex = k;
                            break;
                        }
                    }
                    processingTime[j, m] = processingTimesActivityOrder[j, machineIndex];
                    maxStart += processingTime[j, m];
                }
            }
        }
    }

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

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

        // Interval decisions: time range of jobs on each machine
        // tasks[j][m] is the interval of time of the task of job j 
        // which is processed on machine m
        tasks = new HxExpression[nbJobs, nbMachines];
        for (int j = 0; j < nbJobs; ++j)
        {
            for (int m = 0; m < nbMachines; ++m)
            {
                tasks[j, m] = model.Interval(0, maxStart);

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

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

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

            // Every machine executes a single task at a time
            HxExpression sequenceLambda = model.LambdaFunction(
                i => taskArray[sequence[i], m] < taskArray[sequence[i + 1], m]
            );
            model.Constraint(model.And(model.Range(0, nbJobs - 1), sequenceLambda));
        }

        // Sequence of tasks on each machine
        machinesOrder = new HxExpression[nbJobs];
        for (int j = 0; j < nbJobs; ++j)
        {
            machinesOrder[j] = model.List(nbMachines);
            HxExpression sequence = machinesOrder[j];
            // Every task is scheduled on its corresponding machine
            model.Constraint(model.Count(sequence) == nbMachines);

            // A job has a single task at a time
            HxExpression sequenceLambda = model.LambdaFunction(
                k => taskArray[j, sequence[k]] < taskArray[j, sequence[k + 1]]
            );
            model.Constraint(model.And(model.Range(0, nbMachines - 1), sequenceLambda));
        }

        // Minimize the makespan: end of the last task of the last job
        makespan = model.Max();
        for (int j = 0; j < nbJobs; ++j)
        {
            for (int m = 0; m < nbMachines; ++m)
            {
                makespan.AddOperand(model.End(tasks[j, m]));
            }
        }
        model.Minimize(makespan);

        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))
        {
            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();
            }
        }
        Console.WriteLine("Solution written in file " + fileName);
    }

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: Openshop 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 (Openshop model = new Openshop())
        {
            model.ReadInstance(instanceFile);
            model.Solve(int.Parse(strTimeLimit));
            if (outputFile != null)
                model.WriteSolution(outputFile);
        }
    }
}
Compilation / Execution (Windows)
javac Openshop.java -cp %HX_HOME%\bin\hexaly.jar
java -cp %HX_HOME%\bin\hexaly.jar;. Openshop instances\tai2020_5.txt
Compilation / Execution (Linux)
javac Openshop.java -cp /opt/hexaly_13_0/bin/hexaly.jar
java -cp /opt/hexaly_13_0/bin/hexaly.jar:. Openshop instances/tai2020_5.txt
import java.util.*;
import java.io.*;
import com.hexaly.optimizer.*;

public class Openshop {
    // Number of jobs
    private int nbJobs;
    // Number of machines
    private int nbMachines;
    // Processing time on each machine for each job task
    private long[][] processingTime;
    // 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 : processing order of jobs for each machine
    private HxExpression[] jobsOrder;
    // Decision variables : processing order of machines for each job
    private HxExpression[] machinesOrder;
    // Objective = minimize the makespan: end of the last task of the last job
    private HxExpression makespan;

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

    // The input files follow the "Taillard" format
    public void readInstance(String fileName) throws IOException {
        try (Scanner input = new Scanner(new File(fileName))) {
            input.nextLine();
            nbJobs = input.nextInt();
            nbMachines = input.nextInt();

            input.nextLine();
            input.nextLine();
            // Processing times for each job on each machine
            // (given in the task order, the processing order is a decision variable)
            long[][] processingTimesActivityOrder = new long[nbJobs][nbMachines];
            for (int j = 0; j < nbJobs; ++j) {
                for (int m = 0; m < nbMachines; ++m) {
                    processingTimesActivityOrder[j][m] = input.nextInt();
                }
            }
            // Index of machines for each task
            input.nextLine();
            input.nextLine();
            int[][] machineIndexes = new int[nbJobs][nbMachines];
            for (int j = 0; j < nbJobs; ++j) {
                for (int m = 0; m < nbMachines; ++m) {
                    machineIndexes[j][m] = input.nextInt() - 1;
                }
            }

            // Reorder processing times: processingTime[j][m] is the processing time of the
            // task of job j that is processed on machine m
            processingTime = new long[nbJobs][nbMachines];
            // Trivial upper bound for the start times of the tasks
            maxStart = 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 (machineIndexes[j][k] == m) {
                            machineIndex = k;
                            break;
                        }
                    }
                    processingTime[j][m] = processingTimesActivityOrder[j][machineIndex];
                    maxStart += processingTime[j][m];
                }
            }
        }
    }

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

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

                // Task duration constraints
                model.constraint(model.eq(model.length(tasks[j][m]), processingTime[j][m]));
            }
        }

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

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

            // Every machine executes a single task at a time
            HxExpression mExpr = model.createConstant(m);
            HxExpression sequenceLambda = model
                    .lambdaFunction(i -> model.lt(model.at(taskArray, model.at(sequence, i), mExpr),
                            model.at(taskArray, model.at(sequence, model.sum(i, 1)), mExpr)));
            model.constraint(model.and(model.range(0, nbJobs - 1), sequenceLambda));
        }

        // Sequence of machines for each job
        machinesOrder = new HxExpression[nbJobs];
        for (int j = 0; j < nbJobs; ++j) {
            machinesOrder[j] = model.listVar(nbMachines);
            HxExpression sequence = machinesOrder[j];
            // Every task is scheduled on its corresponding machine
            model.constraint(model.eq(model.count(sequence), nbMachines));

            // A job has a single task at a time
            HxExpression jExpr = model.createConstant(j);
            HxExpression sequenceLambda = model
                    .lambdaFunction(k -> model.lt(model.at(taskArray, jExpr, model.at(sequence, k)),
                            model.at(taskArray, jExpr, model.at(sequence, model.sum(k, 1)))));
            model.constraint(model.and(model.range(0, nbMachines - 1), sequenceLambda));
        }

        // Minimize the makespan: end of the last task
        makespan = model.max();
        for (int m = 0; m < nbMachines; ++m) {
            HxExpression mExpr = model.createConstant(m);
            for (int j = 0; j < nbJobs; ++j) {
                HxExpression jExpr = model.createConstant(j);
                makespan.addOperand(model.end(model.at(taskArray, jExpr, mExpr)));
            }
        }
        model.minimize(makespan);

        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)) {
            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");
            }
        }
        System.out.println("Solution written in file " + fileName);
    }

    public static void main(String[] args) {
        if (args.length < 1) {
            System.out.println("Usage: java Openshop 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()) {
            Openshop model = new Openshop(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);
        }
    }
}