Flexible Job Shop Problem with Machine Changeover Times

import hexaly.optimizer
import sys

# Constant for incompatible machines
INFINITE = 1000000


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

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

    # Number of operations for each job
    nb_operations = [int(lines[j + 1].split()[0]) for j in range(nb_jobs)]

    # Number of tasks
    nb_tasks = sum(nb_operations[j] for j in range(nb_jobs))

    # Processing time for each task, for each machine
    task_processing_time = [[INFINITE for m in range(nb_machines)] for i in range(nb_tasks)]

    # For each job, for each operation, the corresponding task id
    job_operation_task = [[0 for o in range(nb_operations[j])] for j in range(nb_jobs)]

    id = 0
    for j in range(nb_jobs):
        line = lines[j + 1].split()
        tmp = 0
        for o in range(nb_operations[j]):
            nb_machines_operation = int(line[tmp + o + 1])
            for i in range(nb_machines_operation):
                machine = int(line[tmp + o + 2 * i + 2]) - 1
                time = int(line[tmp + o + 2 * i + 3])
                task_processing_time[id][machine] = time
            job_operation_task[j][o] = id
            id = id + 1
            tmp = tmp + 2 * nb_machines_operation

    # Changeover time between two machines
    machine_changeover_time = [[0 for m2 in range(nb_machines)] for m1 in range(nb_machines)]

    for m1 in range(nb_machines):
        line = lines[nb_jobs + 1 + m1].split()
        for m2 in range(nb_machines):
            machine_changeover_time[m1][m2] = int(line[m2])


    # Trivial upper bound for the start times of the tasks
    max_start = sum(
        max(task_processing_time[i][m] for m in range(nb_machines) if task_processing_time[i][m] != INFINITE)
        for i in range(nb_tasks))
    
    return nb_jobs, nb_machines, nb_tasks, task_processing_time, job_operation_task, nb_operations, max_start, machine_changeover_time


def main(instance_file, output_file, time_limit):
    nb_jobs, nb_machines, nb_tasks, task_processing_time_data, job_operation_task, \
        nb_operations, max_start, machine_changeover_time_data = read_instance(instance_file)

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

        # Sequence of tasks on each machine
        jobs_order = [model.list(nb_tasks) for _ in range(nb_machines)]
        machines = model.array(jobs_order)

        # Each task is scheduled on a machine
        model.constraint(model.partition(machines))

        # Only compatible machines can be selected for a task
        for i in range(nb_tasks):
            for m in range(nb_machines):
                if task_processing_time_data[i][m] == INFINITE:
                    model.constraint(model.not_(model.contains(jobs_order[m], i)))

        # For each task, the selected machine
        task_machine = [model.find(machines, i) for i in range(nb_tasks)]

        task_processing_time = model.array(task_processing_time_data)

        # Interval decisions: time range of each task
        tasks = [model.interval(0, max_start) for _ in range(nb_tasks)]

        # The task duration depends on the selected machine
        duration = [model.at(task_processing_time, i, task_machine[i]) for i in range(nb_tasks)]
        for i in range(nb_tasks):
            model.constraint(model.length(tasks[i]) == duration[i])

        task_array = model.array(tasks)

        machine_changeover_time = model.array(machine_changeover_time_data)
        # Precedence constraints between the operations of a job with machine-dependent changeover times
        for j in range(nb_jobs):
            for o in range(nb_operations[j] - 1):
                i1 = job_operation_task[j][o]
                i2 = job_operation_task[j][o + 1]
                model.constraint(model.start(tasks[i2]) >= model.end(tasks[i1]) 
                        + machine_changeover_time[task_machine[i1]][task_machine[i2]])

        # Disjunctive resource constraints between the tasks on a machine
        for m in range(nb_machines):
            sequence = jobs_order[m]
            sequence_lambda = model.lambda_function(
                lambda i: task_array[sequence[i]] < task_array[sequence[i + 1]])
            model.constraint(model.and_(model.range(0, model.count(sequence) - 1), sequence_lambda))

        # Minimize the makespan: end of the last task
        makespan = model.max([model.end(tasks[i]) for i in range(nb_tasks)])
        model.minimize(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 operation of each job, the selected machine, the start and end dates
        if output_file != None:
            with open(output_file, "w") as f:
                print("Solution written in file", output_file)
                for j in range(nb_jobs):
                    for o in range(0, nb_operations[j]):
                        taskIndex = job_operation_task[j][o]
                        f.write(str(j + 1) + "\t" + str(o + 1)
                                + "\t" + str(task_machine[taskIndex].value + 1)
                                + "\t" + str(tasks[taskIndex].value.start())
                                + "\t" + str(tasks[taskIndex].value.end()) + "\n")


if __name__ == '__main__':
    if len(sys.argv) < 2:
        print("Usage: python flexiblejobshop_changeover.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)

#include "optimizer/hexalyoptimizer.h"
#include <algorithm>
#include <fstream>
#include <iostream>
#include <limits>
#include <numeric>
#include <vector>

using namespace hexaly;

class FlexibleJobshop {
private:
    // Number of jobs
    int nbJobs;
    // Number of machines
    int nbMachines;
    // Number of tasks
    int nbTasks;
    // Processing time for each task, for each machine
    std::vector<std::vector<int>> taskProcessingTimeData;
    // Changeover time between two machines
    std::vector<std::vector<int>> machineChangeoverTimeData;
    // For each job, for each operation, the corresponding task id
    std::vector<std::vector<int>> jobOperationTask;
    // Number of operations for each job
    std::vector<int> nbOperations;
    // Trivial upper bound for the start times of the tasks
    int maxStart;
    // Constant for incompatible machines
    const int INFINITE = 1000000;

    // Hexaly Optimizer
    HexalyOptimizer optimizer;
    // Decision variables: time range of each task
    std::vector<HxExpression> tasks;
    // Decision variables: sequence of tasks on each machine
    std::vector<HxExpression> jobsOrder;
    // For each task, the selected machine
    std::vector<HxExpression> taskMachine;
    // Objective = minimize the makespan: end of the last task
    HxExpression makespan;

public:
    FlexibleJobshop() : optimizer() {}

    void readInstance(const std::string& fileName) {
        std::ifstream infile;
        infile.exceptions(std::ifstream::failbit | std::ifstream::badbit);
        infile.open(fileName.c_str());

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

        nbTasks = 0;
        std::vector<std::vector<std::vector<int>>> processingTime = std::vector<std::vector<std::vector<int>>>(nbJobs);
        jobOperationTask.resize(nbJobs);
        nbOperations.resize(nbJobs);
        for (unsigned int j = 0; j < nbJobs; ++j) {
            infile >> nbOperations[j];
            jobOperationTask[j].resize(nbOperations[j]);
            processingTime[j].resize(nbOperations[j]);
            for (unsigned int o = 0; o < nbOperations[j]; ++o) {
                int nbMachinesOperation;
                infile >> nbMachinesOperation;
                taskProcessingTimeData.push_back(std::vector<int>(nbMachines, INFINITE));
                processingTime[j][o].resize(nbMachines, INFINITE);
                for (int m = 0; m < nbMachinesOperation; ++m) {
                    int machine;
                    int time;
                    infile >> machine;
                    infile >> time;
                    processingTime[j][o][machine - 1] = time;
                    taskProcessingTimeData[nbTasks][machine - 1] = time;
                }
                jobOperationTask[j][o] = nbTasks;
                nbTasks += 1;
            }
        }

        machineChangeoverTimeData = std::vector<std::vector<int>>(nbMachines, std::vector<int>(nbMachines));
        for (unsigned int m1 = 0; m1 < nbMachines; ++m1){
            for (unsigned int m2 = 0; m2 < nbMachines; ++m2){
                infile >> machineChangeoverTimeData[m1][m2];
            }
        }
        infile.close();

        // Trivial upper bound for the start times of the tasks
        maxStart = 0;
        for (unsigned int j = 0; j < nbJobs; ++j) {
            for (unsigned int o = 0; o < nbOperations[j]; ++o) {
                int maxProcessingTime = 0;
                for (unsigned int m = 0; m < nbMachines; ++m) {
                    if (processingTime[j][o][m] != INFINITE && processingTime[j][o][m] > maxProcessingTime)
                        maxProcessingTime = processingTime[j][o][m];
                }
                maxStart += maxProcessingTime;
            }
        }
    }

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

        // Sequence of tasks on each machine
        jobsOrder.resize(nbMachines);
        HxExpression machines = model.array();
        for (unsigned int m = 0; m < nbMachines; ++m) {
            jobsOrder[m] = model.listVar(nbTasks);
            machines.addOperand(jobsOrder[m]);
        }

        // Each task is scheduled on a machine
        model.constraint(model.partition(machines));

        // Only compatible machines can be selected for a task
        for (int i = 0; i < nbTasks; ++i) {
            for (unsigned int m = 0; m < nbMachines; ++m) {
                if (taskProcessingTimeData[i][m] == INFINITE) {
                    model.constraint(!model.contains(jobsOrder[m], i));
                }
            }
        }

        taskMachine.resize(nbTasks);
        HxExpression taskProcessingTime = model.array();
        for (int i = 0; i < nbTasks; ++i) {
            // For each task, the selected machine
            taskMachine[i] = model.find(machines, i);
            taskProcessingTime.addOperand(
                model.array(taskProcessingTimeData[i].begin(), taskProcessingTimeData[i].end()));
        }

        tasks.resize(nbTasks);
        std::vector<HxExpression> duration(nbTasks);
        for (int i = 0; i < nbTasks; ++i) {
            // Interval decisions: time range of each task
            tasks[i] = model.intervalVar(0, maxStart);

            // The task duration depends on the selected machine
            duration[i] = model.at(taskProcessingTime, i, taskMachine[i]);
            model.constraint(model.length(tasks[i]) == duration[i]);
        }
        HxExpression taskArray = model.array(tasks.begin(), tasks.end());


        HxExpression machineChangeoverTime = model.array();
        for (int m1 = 0; m1 < nbMachines; ++m1) {
            machineChangeoverTime.addOperand(
                    model.array(machineChangeoverTimeData[m1].begin(), machineChangeoverTimeData[m1].end()));
            }

        // Precedence constraints between the operations of a job with machine-dependent changeover times
        for (unsigned int j = 0; j < nbJobs; ++j) {
            for (unsigned int o = 0; o < nbOperations[j] - 1; ++o) {
                int i1 = jobOperationTask[j][o];
                int i2 = jobOperationTask[j][o + 1];
                model.constraint(model.start(tasks[i2]) >= model.end(tasks[i1])
                        + machineChangeoverTime[taskMachine[i1]][taskMachine[i2]]);
            }
        }

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

        // Minimize the makespan: end of the last task
        makespan = model.max();
        for (int i = 0; i < nbTasks; ++i) {
            makespan.addOperand(model.end(tasks[i]));
        }
        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 operation of each job, the selected machine, the start and end dates */
    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);
        }
        std::cout << "Solution written in file " << fileName << std::endl;

        for (unsigned int j = 0; j < nbJobs; ++j) {
            for (unsigned int o = 0; o < nbOperations[j]; ++o) {
                int taskIndex = jobOperationTask[j][o];
                outfile << j + 1 << "\t" << o + 1 << "\t" << taskMachine[taskIndex].getValue() + 1 << "\t"
                        << tasks[taskIndex].getIntervalValue().getStart() << "\t"
                        << tasks[taskIndex].getIntervalValue().getEnd() << std::endl;
            }
        }
        outfile.close();
    }
};

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

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

    FlexibleJobshop model;
    try {
        model.readInstance(instanceFile);
        const int timeLimit = atoi(strTimeLimit);
        model.solve(timeLimit);
        if (outputFile != NULL)
            model.writeSolution(outputFile);
        return 0;
    } catch (const std::exception& e) {
        std::cerr << "An error occurred: " << e.what() << std::endl;
        return 1;
    }
}

using System;
using System.IO;
using System.Linq;
using Hexaly.Optimizer;

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

    // Number of machines
    private int nbMachines;

    // Number of tasks
    private int nbTasks;

    // Processing time for each task, for each machine
    private long[][] taskProcessingTimeData;

    // Changeover time between two machines
    private int[][] machineChangeoverTimeData;

    // For each job, for each operation, the corresponding task id
    private int[][] jobOperationTask;

    // Number of operations for each job;
    private int[] nbOperations;

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

    // Constant for incompatible machines
    private const long INFINITE = 1000000;

    // 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;

    // For each task, the selected machine
    private HxExpression[] taskMachine;

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

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

    public void ReadInstance(string fileName)
    {
        using (StreamReader input = new StreamReader(fileName))
        {
            char[] separators = new char[] { '\t', ' ' };
            string[] splitted = input
                .ReadLine()
                .Split(separators, StringSplitOptions.RemoveEmptyEntries);
            nbJobs = int.Parse(splitted[0]);
            nbMachines = int.Parse(splitted[1]);

            nbTasks = 0;
            long[][][] processingTime = new long[nbJobs][][];
            jobOperationTask = new int[nbJobs][];
            nbOperations = new int[nbJobs];
            for (int j = 0; j < nbJobs; ++j)
            {
                splitted = input
                    .ReadLine()
                    .Split(separators, StringSplitOptions.RemoveEmptyEntries);
                nbOperations[j] = int.Parse(splitted[0]);
                jobOperationTask[j] = new int[nbOperations[j]];
                processingTime[j] = new long[nbOperations[j]][];
                int k = 1;
                for (int o = 0; o < nbOperations[j]; ++o)
                {
                    int nbMachinesOperation = int.Parse(splitted[k]);
                    k++;
                    processingTime[j][o] = Enumerable.Repeat((long)INFINITE, nbMachines).ToArray();
                    for (int m = 0; m < nbMachinesOperation; ++m)
                    {
                        int machine = int.Parse(splitted[k]) - 1;
                        long time = long.Parse(splitted[k + 1]);
                        processingTime[j][o][machine] = time;
                        k += 2;
                    }
                    jobOperationTask[j][o] = nbTasks;
                    nbTasks++;
                }
            }

            machineChangeoverTimeData = new int[nbMachines][];
            for (int m1 = 0; m1 < nbMachines; ++m1)
            {
                machineChangeoverTimeData[m1] = new int[nbTasks];
                    splitted = input.
                        ReadLine()
                        .Split(separators, StringSplitOptions.RemoveEmptyEntries);

                machineChangeoverTimeData[m1] = new int[nbMachines];
                for (int m2 = 0; m2 < nbMachines; ++m2)
                {
                    machineChangeoverTimeData[m1][m2] = int.Parse(splitted[m2]);
                }
            }

            // Trivial upper bound for the start times of the tasks
            maxStart = 0;
            taskProcessingTimeData = new long[nbTasks][];
            for (int j = 0; j < nbJobs; ++j)
            {
                long maxProcessingTime = 0;
                for (int o = 0; o < nbOperations[j]; ++o)
                {
                    int task = jobOperationTask[j][o];
                    taskProcessingTimeData[task] = new long[nbMachines];
                    for (int m = 0; m < nbMachines; ++m)
                    {
                        taskProcessingTimeData[task][m] = processingTime[j][o][m];
                        if (
                            processingTime[j][o][m] != INFINITE
                            && processingTime[j][o][m] > maxProcessingTime
                        )
                        {
                            maxProcessingTime = processingTime[j][o][m];
                        }
                    }
                    maxStart += maxProcessingTime;
                }
            }
        }
    }

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

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

        // Sequence of tasks on each machine
        jobsOrder = new HxExpression[nbMachines];
        HxExpression machines = model.Array();
        for (int m = 0; m < nbMachines; ++m)
        {
            jobsOrder[m] = model.List(nbTasks);
            machines.AddOperand(jobsOrder[m]);
        }

        // Each task is scheduled on a machine
        model.Constraint(model.Partition(machines));

        // Only compatible machines can be selected for a task
        for (int i = 0; i < nbTasks; ++i)
        {
            for (int m = 0; m < nbMachines; ++m)
            {
                if (taskProcessingTimeData[i][m] == INFINITE)
                    model.Constraint(!model.Contains(jobsOrder[m], i));
            }
        }

        // For each task, the selected machine
        taskMachine = new HxExpression[nbTasks];
        for (int i = 0; i < nbTasks; ++i)
        {
            taskMachine[i] = model.Find(machines, i);
        }

        tasks = new HxExpression[nbTasks];
        HxExpression[] duration = new HxExpression[nbTasks];
        HxExpression taskProcessingTime = model.Array(taskProcessingTimeData);
        for (int i = 0; i < nbTasks; ++i)
        {
            // Interval decisions: time range of each task
            tasks[i] = model.Interval(0, maxStart);

            // The task duration depends on the selected machine
            HxExpression iExpr = model.CreateConstant(i);
            duration[i] = model.At(taskProcessingTime, iExpr, taskMachine[i]);
            model.Constraint(model.Length(tasks[i]) == duration[i]);
        }

        HxExpression taskArray = model.Array(tasks);
        HxExpression machineChangeoverTime = model.Array(machineChangeoverTimeData);
        // Precedence constraints between the operations of a job with machine-dependent changeover times
        for (int j = 0; j < nbJobs; ++j)
        {
            for (int o = 0; o < nbOperations[j] - 1; ++o)
            {
                int i1 = jobOperationTask[j][o];
                int i2 = jobOperationTask[j][o + 1];
                model.Constraint(model.Start(tasks[i2]) >= model.End(tasks[i1])
                        + machineChangeoverTime[taskMachine[i1]][taskMachine[i2]]); 
            }
        }

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

        // Minimize the makespan: end of the last task
        makespan = model.Max();
        for (int i = 0; i < nbTasks; ++i)
        {
            makespan.AddOperand(model.End(tasks[i]));
        }
        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 operation of each job, the selected machine, the start and end dates */
    public void WriteSolution(string fileName)
    {
        using (StreamWriter output = new StreamWriter(fileName))
        {
            Console.WriteLine("Solution written in file " + fileName);
            for (int j = 1; j <= nbJobs; ++j)
            {
                for (int o = 1; o <= nbOperations[j - 1]; ++o)
                {
                    int taskIndex = jobOperationTask[j - 1][o - 1];
                    output.WriteLine(
                        j
                            + "\t"
                            + o
                            + "\t"
                            + taskMachine[taskIndex].GetValue()
                            + "\t"
                            + tasks[taskIndex].GetIntervalValue().Start()
                            + "\t"
                            + tasks[taskIndex].GetIntervalValue().End()
                    );
                }
            }
        }
    }

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: FlexibleJobshopChangeover 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 (FlexibleJobshopChangeover model = new FlexibleJobshopChangeover())
        {
            model.ReadInstance(instanceFile);
            model.Solve(int.Parse(strTimeLimit));
            if (outputFile != null)
                model.WriteSolution(outputFile);
        }
    }
}

import java.util.*;
import java.io.*;
import com.hexaly.optimizer.*;

public class FlexibleJobshopChangeover {
    // Number of jobs
    private int nbJobs;
    // Number of machines
    private int nbMachines;
    // Number of tasks
    private int nbTasks;
    // Processing time for each task, for each machine
    private long[][] taskProcessingTimeData;
    // Changeover time between two machines
    private int[][] machineChangeoverTimeData;
    // For each job, for each operation, the corresponding task id
    private int[][] jobOperationTask;
    // Number of operations for each job;
    private int[] nbOperations;
    // Trivial upper bound for the start times of the tasks
    private long maxStart;
    // Constant for incompatible machines
    private final int INFINITE = 1000000;

    // 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;
    // For each task, the selected machine
    private HxExpression[] taskMachine;
    // Objective = minimize the makespan: end of the last task
    private HxExpression makespan;

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

    public void readInstance(String fileName) throws IOException {
        try (Scanner input = new Scanner(new File(fileName))) {
            nbJobs = input.nextInt();
            nbMachines = input.nextInt();
            input.next(); // skip last number

            nbTasks = 0;
            long[][][] processingTime = new long[nbJobs][][];
            jobOperationTask = new int[nbJobs][];
            nbOperations = new int[nbJobs];
            for (int j = 0; j < nbJobs; ++j) {
                nbOperations[j] = input.nextInt();
                jobOperationTask[j] = new int[nbOperations[j]];
                processingTime[j] = new long[nbOperations[j]][nbMachines];
                for (int o = 0; o < nbOperations[j]; ++o) {
                    int nbMachinesOperation = input.nextInt();
                    Arrays.fill(processingTime[j][o], INFINITE);
                    for (int m = 0; m < nbMachinesOperation; ++m) {
                        int machine = input.nextInt() - 1;
                        long time = input.nextLong();
                        processingTime[j][o][machine] = time;
                    }
                    jobOperationTask[j][o] = nbTasks;
                    nbTasks++;
                }
            }


            machineChangeoverTimeData = new int[nbMachines][nbMachines];
            for (int m1 = 0; m1 < nbMachines; ++m1) {
                for (int m2 = 0; m2 < nbMachines; ++m2) {
                    machineChangeoverTimeData[m1][m2] = input.nextInt();
                }
            }

            // Trivial upper bound for the start times of the tasks
            maxStart = 0;
            taskProcessingTimeData = new long[nbTasks][];
            for (int j = 0; j < nbJobs; ++j) {
                long maxProcessingTime = 0;
                for (int o = 0; o < nbOperations[j]; ++o) {
                    int task = jobOperationTask[j][o];
                    taskProcessingTimeData[task] = new long[nbMachines];
                    for (int m = 0; m < nbMachines; ++m) {
                        taskProcessingTimeData[task][m] = processingTime[j][o][m];
                        if (processingTime[j][o][m] != INFINITE && processingTime[j][o][m] > maxProcessingTime) {
                            maxProcessingTime = processingTime[j][o][m];
                        }
                    }
                    maxStart += maxProcessingTime;
                }
            }
        }
    }

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

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

        // Each task is scheduled on a machine
        model.constraint(model.partition(machines));

        // Only compatible machines can be selected for a task
        for (int i = 0; i < nbTasks; ++i) {
            for (int m = 0; m < nbMachines; ++m) {
                if (taskProcessingTimeData[i][m] == INFINITE) {
                    model.constraint(model.not(model.contains(jobsOrder[m], i)));
                }
            }
        }

        // For each task, the selected machine
        taskMachine = new HxExpression[nbTasks];
        for (int i = 0; i < nbTasks; ++i) {
            taskMachine[i] = model.find(machines, i);
        }

        HxExpression taskProcessingTime = model.array(taskProcessingTimeData);

        tasks = new HxExpression[nbTasks];
        HxExpression[] duration = new HxExpression[nbTasks];
        for (int i = 0; i < nbTasks; ++i) {
            // Interval decisions: time range of each task
            tasks[i] = model.intervalVar(0, maxStart);

            // The task duration depends on the selected machine
            HxExpression iExpr = model.createConstant(i);
            duration[i] = model.at(taskProcessingTime, iExpr, taskMachine[i]);
            model.constraint(model.eq(model.length(tasks[i]), duration[i]));
        }

        HxExpression taskArray = model.array(tasks);
        HxExpression machineChangeoverTime = model.array(machineChangeoverTimeData);
    
        // Precedence constraints between the operations of a job with machine-dependent changeover times
        for (int j = 0; j < nbJobs; ++j) {
            for (int o = 0; o < nbOperations[j] - 1; ++o) {
                int i1 = jobOperationTask[j][o];
                int i2 = jobOperationTask[j][o + 1];
                model.constraint(model.geq(model.start(tasks[i2]), model.sum(model.end(tasks[i1]), 
                        model.at(machineChangeoverTime, taskMachine[i1], taskMachine[i2]))));
            }
        }

        // Disjunctive resource constraints between the tasks on a machine
        for (int m = 0; m < nbMachines; ++m) {
            HxExpression sequence = jobsOrder[m];
            HxExpression sequenceLambda = model.lambdaFunction(i -> model
                    .lt(model.at(taskArray, model.at(sequence, i)),
                            model.at(taskArray, model.at(sequence, model.sum(i, 1)))));
            model.constraint(model.and(model.range(0, model.sub(model.count(sequence), 1)), sequenceLambda));
        }

        // Minimize the makespan: end of the last task
        makespan = model.max();
        for (int i = 0; i < nbTasks; ++i) {
            makespan.addOperand(model.end(tasks[i]));
        }
        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 operation of each job, the selected machine, the start and end
     * dates
     */
    public void writeSolution(String fileName) throws IOException {
        try (PrintWriter output = new PrintWriter(fileName)) {
            System.out.println("Solution written in file " + fileName);

            for (int j = 1; j <= nbJobs; ++j) {
                for (int o = 1; o <= nbOperations[j - 1]; ++o) {
                    int taskIndex = jobOperationTask[j - 1][o - 1];
                    output.write(j + "\t" + o
                            + "\t" + taskMachine[taskIndex].getValue()
                            + "\t" + tasks[taskIndex].getIntervalValue().getStart()
                            + "\t" + tasks[taskIndex].getIntervalValue().getEnd() + "\n");
                }
            }
        }
    }

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