Job Shop Scheduling Problem with Intensity

import hexaly.optimizer
import sys


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

    first_line = lines[1].split()
    # Number of jobs
    nb_jobs = int(first_line[0])
    # Number of machines
    nb_machines = int(first_line[1])
    # Time horizon: number of time steps
    time_horizon = int(first_line[2])

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

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

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

    # Intensity for each machine for each time step
    intensity = [[int(lines[i].split()[j]) for j in range(time_horizon)]
                 for i in range(5 + 2 * nb_jobs, 5 + 2 * nb_jobs + nb_machines)]

    return nb_jobs, nb_machines, time_horizon, processing_time, machine_order, intensity


def main(instance_file, output_file, time_limit):
    nb_jobs, nb_machines, time_horizon, processing_time, \
        machine_order, intensity = 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, time_horizon) for m in range(nb_machines)]
                 for j in range(nb_jobs)]

        # Create Hexaly arrays to be able to access them with "at" operators
        task_array = model.array(tasks)
        intensity_array = model.array(intensity)

        # The sum of the machine's intensity over the duration of the task must be
        # greater than its processing time
        for m in range(nb_machines):
            intensity_lambda = model.lambda_function(lambda t:
                                                     model.at(intensity_array, m, t))
            for j in range(nb_jobs):
                model.constraint(
                    model.sum(tasks[j][m], intensity_lambda)
                    >= processing_time[j][m])

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

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

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

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

        # Minimize the makespan: end of the last task of the last job
        makespan = model.max([model.end(tasks[j][machine_order[j][nb_machines - 1]])
                              for j in range(nb_jobs)])
        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 machine, the job sequence
        #
        if output_file != None:
            final_jobs_order = [list(jobs_order[m].value) for m in range(nb_machines)]
            with open(output_file, "w") as f:
                print("Solution written in file ", output_file)
                for m in range(nb_machines):
                    for j in range(nb_jobs):
                        f.write(str(final_jobs_order[m][j]) + " ")
                    f.write("\n")


if __name__ == '__main__':
    if len(sys.argv) < 2:
        print("Usage: python jobshop_intensity.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 10
    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 JobshopIntensity {
private:
    // Number of jobs
    int nbJobs;
    // Number of machines
    int nbMachines;
    // Number of time steps
    int timeHorizon;

    // Processing order of machines for each job
    std::vector<std::vector<int>> machineOrder;
    // Processing time on each machine for each job (given in the machine order)
    std::vector<std::vector<int>> processingTime;
    // Intensity of each machine for each time step
    std::vector<std::vector<int>> intensity;

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

public:
    JobshopIntensity() : optimizer() {}

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

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

        // Processing times for each job on each machine (given in the processing order)
        infile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
        std::vector<std::vector<int>> processingTimesInProcessingOrder =
            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 >> processingTimesInProcessingOrder[j][m];
            }
        }

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

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

        // Intensity for each machine for each time step
        infile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
        infile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
        intensity.resize(nbMachines);
        for (int m = 0; m < nbMachines; ++m) {
            intensity[m].resize(timeHorizon);
            for (int t = 0; t < timeHorizon; ++t) {
                int x;
                infile >> x;
                intensity[m][t] = x;
            }
        }

        infile.close();
    }

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

        // 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.resize(nbJobs);
        for (unsigned int j = 0; j < nbJobs; ++j) {
            tasks[j].resize(nbMachines);
            for (unsigned int m = 0; m < nbMachines; ++m) {
                tasks[j][m] = model.intervalVar(0, timeHorizon);
            }
        }

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

        HxExpression intensityArray = model.array();
        for (int m = 0; m < nbMachines; ++m) {
            intensityArray.addOperand(model.array(intensity[m].begin(), intensity[m].end()));
        }

        // The sum of the machine's intensity over the duration of the task must be greater than its processing time
        for (int m = 0; m < nbMachines; ++m) {
            HxExpression intensityLambda =
                model.createLambdaFunction([&](HxExpression t) { return model.at(intensityArray, m, t); });
            for (int j = 0; j < nbJobs; ++j) {
                model.constraint(model.sum(model.at(taskArray, j, m), intensityLambda) >= processingTime[j][m]);
            }
        }

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

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

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

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

        // Minimize the makespan: end of the last task of the last job
        makespan = model.max();
        for (int j = 0; j < nbJobs; ++j) {
            makespan.addOperand(model.end(tasks[j][machineOrder[j][nbMachines - 1]]));
        }
        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 */
    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 (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();
    }
};

int main(int argc, char** argv) {
    if (argc < 2) {
        std::cout << "Usage: jobshop_intensity 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] : "10";

    JobshopIntensity 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 Hexaly.Optimizer;

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

    // Number of machines
    private int nbMachines;

    // Number of time steps
    private long timeHorizon;

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

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

    // Intensity of each machine for each time step
    private int[,] intensity;

    // Hexaly Optimizer
    private HexalyOptimizer optimizer;

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

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

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

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

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

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

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

            // Reorder processing times: processingTime[j, m] is the processing time of the
            // task of job j that is processed on machine m
            processingTime = new long[nbJobs, nbMachines];
            for (int j = 0; j < nbJobs; ++j)
            {
                for (int m = 0; m < nbMachines; ++m)
                {
                    int machineIndex = nbMachines;
                    for (int k = 0; k < nbMachines; ++k)
                    {
                        if (machineOrder[j, k] == m)
                        {
                            machineIndex = k;
                            break;
                        }
                    }
                    processingTime[j, m] = processingTimesInProcessingOrder[j, machineIndex];
                }
            }

            // Intensity for each machine for each time step
            input.ReadLine();
            intensity = new int[nbMachines, timeHorizon];
            for (int m = 0; m < nbMachines; ++m)
            {
                splitted = input.ReadLine().Trim().Split(' ');
                for (int t = 0; t < timeHorizon; ++t)
                    intensity[m, t] = int.Parse(splitted[t]);
            }
        }
    }

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

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

        // Interval decisions: time range of each task
        // tasks[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, timeHorizon);
            }
        }

        // Create HexalyOptimizer arrays to be able to access them with "at" operators
        HxExpression taskArray = model.Array(tasks);
        HxExpression intensityArray = model.Array(intensity);

        // The sum of the machine's intensity over the duration of the task must be greater than its processing time
        for (int m = 0; m < nbMachines; ++m)
        {
            HxExpression intensityLambda = model.LambdaFunction(t => intensityArray[m, t]);
            for (int j = 0; j < nbJobs; ++j)
                model.Constraint(model.Sum(taskArray[j, m], intensityLambda)
                        >= processingTime[j, m]);
        }

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

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

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

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

        // Minimize the makespan: end of the last task of the last job
        makespan = model.Max();
        for (int j = 0; j < nbJobs; ++j)
            makespan.AddOperand(model.End(tasks[j, machineOrder[j, nbMachines - 1]]));
        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))
        {
            Console.WriteLine("Solution written in file " + fileName);
            for (int m = 0; m < nbMachines; ++m)
            {
                HxCollection finalJobsOrder = jobsOrder[m].GetCollectionValue();
                for (int i = 0; i < nbJobs; i++)
                {
                    int j = (int)finalJobsOrder.Get(i);
                    output.Write(j + " ");
                }
                output.WriteLine();
            }
        }
    }

    public static void Main(string[] args)
    {
        if (args.Length < 1)
        {
            Console.WriteLine("Usage: JobshopIntensity 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] : "10";

        using (JobshopIntensity model = new JobshopIntensity())
        {
            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 JobshopIntensity {
    // Number of jobs
    private int nbJobs;
    // Number of machines
    private int nbMachines;
    // Time horizon: number of time steps
    private int timeHorizon;
    // Processing time on each machine for each job (given in the machine order)
    private long[][] processingTime;
    // Processing order of machines for each job
    private int[][] machineOrder;
    // Intensity of each machine for each time step
    private int[][] intensity;

    // Hexaly Optimizer
    final HexalyOptimizer optimizer;
    // Decision variables: time range of each task
    private HxExpression[][] tasks;
    // Decision variables: sequence of tasks on each machine
    private HxExpression[] jobsOrder;
    // Objective = minimize the makespan: end of the last task of the last job
    private HxExpression makespan;

    public JobshopIntensity(HexalyOptimizer optimizer, String fileName) throws IOException {
        this.optimizer = optimizer;
    }

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

            input.nextLine();
            input.nextLine();
            // Processing times for each job on each machine (given in the processing order)
            long[][] processingTimesInProcessingOrder = new long[nbJobs][nbMachines];
            for (int j = 0; j < nbJobs; ++j) {
                for (int m = 0; m < nbMachines; ++m) {
                    processingTimesInProcessingOrder[j][m] = input.nextInt();
                }
            }
            // Processing order of machines for each job
            input.nextLine();
            input.nextLine();
            machineOrder = new int[nbJobs][nbMachines];
            for (int j = 0; j < nbJobs; ++j) {
                for (int m = 0; m < nbMachines; ++m) {
                    machineOrder[j][m] = input.nextInt() - 1;
                }
            }
            // Reorder processing times: processingTime[j][m] is the processing time of the
            // task of job j that is processed on machine m
            processingTime = new long[nbJobs][nbMachines];
            for (int j = 0; j < nbJobs; ++j) {
                for (int m = 0; m < nbMachines; ++m) {
                    int machineIndex = nbMachines;
                    for (int k = 0; k < nbMachines; ++k) {
                        if (machineOrder[j][k] == m) {
                            machineIndex = k;
                            break;
                        }
                    }
                    processingTime[j][m] = processingTimesInProcessingOrder[j][machineIndex];
                }
            }

            // Intensity for each machine for each time step
            input.nextLine();
            input.nextLine();
            intensity = new int[nbMachines][timeHorizon];
            for (int m = 0; m < nbMachines; ++m) {
                for (int t = 0; t < timeHorizon; ++t) {
                    intensity[m][t] = input.nextInt();
                }
            }
        }
    }

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

        // 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 = new HxExpression[nbJobs][nbMachines];
        for (int j = 0; j < nbJobs; ++j) {
            for (int m = 0; m < nbMachines; ++m) {
                tasks[j][m] = model.intervalVar(0, timeHorizon);
            }
        }

        // Create HexalyOptimizer arrays to be able to access them with "at" operators
        HxExpression taskArray = model.array(tasks);
        HxExpression intensityArray = model.array(intensity);

        // The sum of the machine's intensity over the duration of the task must be
        // greater than its processing time
        for (int m = 0; m < nbMachines; ++m) {
            HxExpression mExpr = model.createConstant(m);
            HxExpression intensityLambda = model.lambdaFunction(
                    t -> model.at(intensityArray, mExpr, t));
            for (int j = 0; j < nbJobs; ++j) {
                HxExpression jExpr = model.createConstant(j);
                model.constraint(model.geq(model.sum(
                        model.at(taskArray, jExpr, mExpr),
                        intensityLambda), processingTime[j][m]));
            }
        }

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

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

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

            // Disjunctive resource constraints between the tasks on a machine
            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));
        }

        // Minimize the makespan: end of the last task of the last job
        makespan = model.max();
        for (int j = 0; j < nbJobs; ++j) {
            makespan.addOperand(model.end(tasks[j][machineOrder[j][nbMachines - 1]]));
        }
        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)) {
            System.out.println("Solution written in file " + fileName);

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

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

        try (HexalyOptimizer optimizer = new HexalyOptimizer()) {
            JobshopIntensity model = new JobshopIntensity(optimizer, instanceFile);
            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);
        }
    }
}