Open Shop Scheduling Problem (OSSP)

Problem

In the Open Shop Scheduling Problem (OSSP), a set of jobs has to be processed on every machine in the shop. Each job consists in an unordered sequence 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 cannot start an activity while another activity of the job is still running. 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.

Principles learned

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

The Hexaly model for the Open Shop Scheduling Problem (OSSP) uses interval decision variables to represent the time ranges of the activities. The length of the interval is constrained by each activity’s processing time.

In addition to intervals, we also use list decision variables. As in the Job Shop example, a list models the ordering of activities on a machine or within a job. Using the ‘count’ operator to constrain the size of the lists, we ensure that each job is processed on each machine.

The disjunctive resource constraints — each machine can only process one activity at a time — can be reformulated 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.

We use the same strategy to model the disjunctive activity constraints. For all jobs and all i, the activity in position i+1 for a job must start after the end of the activity in position i for this job. Similarly to the disjunctive resource constraints, we model these constraints with a lambda function used within a variadic ‘and’ operator over all activities constituting each job.

The objective consists in minimizing the makespan, which is the time when all the activities have been processed.

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);
        }
    }
}