Job Shop Scheduling Problem with Intensity

Problem

In the Job Shop Scheduling Problem with Intensity, a set of jobs has to be processed on every machine in the shop. Each job consists in an ordered sequence of tasks (called activities), each representing the processing of the job on one of the machines. The jobs have one activity per machine, and cannot start an activity until the previous activity of the job has been completed. Each activity has a given processing time, and each machine can only process one activity at a time. The intensity with which machines process activities varies over time. In particular, an intensity of 0 means that the machine is off. At each time step, each ongoing activity’s progress increases by the current intensity of its machine. An activity starting at time t is considered completed once the sum of its machine’s intensity since time t reaches its processing 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

We provide instances from the FT and LA datasets for the Job Shop Scheduling Problem, with additional intensity matrices. The format of the data files is as follows:

  • First line: number of jobs, number of machines, and time horizon
  • For each job: the processing time on each machine (given in the processing order)
  • For each job: the processing order (ordered list of visited machines)
  • For each machine: its intensity for each time step.

Program

The Hexaly model for the Job Shop Scheduling Problem with Intensity is very similar to the model for the Job Shop Scheduling Problem. The original decision variables remain unchanged. On the one hand, we use interval decision variables to model the time ranges of the activities. On the other hand, we use list decision variables, representing the order of the activities scheduled on the machines.

The intensity constraints define the relationship between the time ranges of the activities and the intensity of the machines over time: the sum of intensities over the span of an activity must be greater than its processing time. To model these constraints, we use a lambda function within a ‘sum’ operator.

The precedence and disjunctive resource constraints are modeled in the same way as for the Job Shop Scheduling Problem. The makespan to be minimized is the time when all jobs have been processed.

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


function input() {
    local usage = "Usage: hexaly jobshop_intensity.hxm inFileName=instanceFile "
            + "[outFileName=outputFile] [hxTimeLimit=timeLimit]";
    if (inFileName == nil) throw usage;

    inFile = io.openRead(inFileName);
    inFile.readln();
    nbJobs = inFile.readInt();
    nbMachines = inFile.readInt();
    timeHorizon = inFile.readInt();
    inFile.readln();

    // Processing times for each job on each machine (given in the processing order)
    processingTimesInProcessingOrder[j in 0...nbJobs][m in 0...nbMachines] = inFile.readInt();
    inFile.readln();
    for [j in 0...nbJobs][k in 0...nbMachines] {
        local m = inFile.readInt() - 1;
        // Processing order of machines for each job
        machineOrder[j][k] = m;
        // Reorder processing times: processingTime[j][m] is the processing time of the
        // task of job j that is processed on machine m
        processingTime[j][m] = processingTimesInProcessingOrder[j][k];
    }

    inFile.readln();
    // Intensity for each machine for each time step
    intensity[m in 0...nbMachines][t in 0...timeHorizon] = inFile.readInt();

    inFile.close();
}

function 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[j in 0...nbJobs][m in 0...nbMachines] <- interval(0, timeHorizon);

    // The sum of the machine's intensity over the duration of the task must be
    // greater than its processing time
    for [j in 0...nbJobs][m in 0...nbMachines] {
        constraint sum(tasks[j][m], t => intensity[m][t]) >= processingTime[j][m];
    }

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

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

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

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

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

    minimize makespan;
}

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

/* Write the solution in a file with the following format:
 *  - for each machine, the job sequence */
function output() {
    if (outFileName != nil) {
        outFile = io.openWrite(outFileName);
        println("Solution written in file ", outFileName);
        for [m in 0...nbMachines]
            outFile.println[j in 0...nbJobs](jobsOrder[m].value[j], " ");
    }
}
Execution (Windows)
set PYTHONPATH=%HX_HOME%\bin\python
python jobshop_intensity.py instances\i01_ft06.txt
Execution (Linux)
export PYTHONPATH=/opt/hexaly_13_0/bin/python
python jobshop_intensity.py instances/i01_ft06.txt
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)
Compilation / Execution (Windows)
cl /EHsc jobshop_intensity.cpp -I%HX_HOME%\include /link %HX_HOME%\bin\hexaly130.lib
jobshop_intensity instances\i01_ft06.txt
Compilation / Execution (Linux)
g++ jobshop_intensity.cpp -I/opt/hexaly_13_0/include -lhexaly130 -lpthread -o jobshop_intensity
./jobshop_intensity instances/i01_ft06.txt
#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;
    }
}
Compilation / Execution (Windows)
copy %HX_HOME%\bin\Hexaly.NET.dll .
csc JobshopIntensity.cs /reference:Hexaly.NET.dll
JobshopIntensity instances\i01_ft06.txt
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);
        }
    }
}
Compilation / Execution (Windows)
javac JobshopIntensity.java -cp %HX_HOME%\bin\hexaly.jar
java -cp %HX_HOME%\bin\hexaly.jar;. JobshopIntensity instances\i01_ft06.txt
Compilation / Execution (Linux)
javac JobshopIntensity.java -cp /opt/hexaly_13_0/bin/hexaly.jar
java -cp /opt/hexaly_13_0/bin/hexaly.jar:. JobshopIntensity instances/i01_ft06.txt
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);
        }
    }
}