Job Shop Scheduling Problem with Intensity¶

Principles learned¶

  • Add multiple list decision variables

  • Constrain the number of elements in a list

  • Use interval decision variables

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

Problem¶

../_images/jobshop_intensity.svg

In the job shop scheduling problem with intensity, a set of jobs has to be processed on every machine of 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. Each job has one activity per machine, and cannot start an activity while the previous activity of the job is not 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.

Download the example


Data¶

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 model is very similar to the original Job Shop Problem, to which we add intensity constraints. The original decision variables remain unchanged: interval decision variables to model the time ranges of the activities, and a list decision variable for each machine, representing the order of the activities scheduled on this machine.

We start by writing the intensity constraints, which define the relationship between the time ranges of the activities and the intensity of the machines over time: the sum of intensities over the duration 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 original Job Shop problem. The makespan to be minimized is the time when all jobs have been processed.

Execution:
localsolver jobshop_intensity.lsp inFileName=instances/i01_ft06.txt [outFileName=] [lsTimeLimit=]
use io;


function input() {
    local usage = "Usage: localsolver jobshop_intensity.lsp inFileName=instanceFile "
            + "[outFileName=outputFile] [lsTimeLimit=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 (lsTimeLimit == nil) lsTimeLimit = 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=%LS_HOME%\bin\python
python jobshop_intensity.py instances\i01_ft06.txt
Execution (Linux)
export PYTHONPATH=/opt/localsolver_13_0/bin/python
python jobshop_intensity.py instances/i01_ft06.txt
import localsolver
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 localsolver.LocalSolver() as ls:
        #
        # Declare the optimization model
        #
        model = ls.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 LocalSolver 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 solver
        ls.param.time_limit = time_limit

        ls.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%LS_HOME%\include /link %LS_HOME%\bin\localsolver130.lib
jobshop_intensity instances\i01_ft06.txt
Compilation / Execution (Linux)
g++ jobshop_intensity.cpp -I/opt/localsolver_13_0/include -llocalsolver130 -lpthread -o jobshop_intensity
./jobshop_intensity instances/i01_ft06.txt
#include "localsolver.h"
#include <algorithm>
#include <fstream>
#include <iostream>
#include <limits>
#include <numeric>
#include <vector>

using namespace localsolver;

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;

    // LocalSolver
    LocalSolver localsolver;
    // Decision variables: time range of each task
    std::vector<std::vector<LSExpression>> tasks;
    // Decision variables: sequence of tasks on each machine
    std::vector<LSExpression> jobsOrder;
    // Objective = minimize the makespan: end of the last task of the last job
    LSExpression makespan;

public:
    JobshopIntensity() : localsolver() {}

    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
        LSModel model = localsolver.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 LocalSolver arrays to be able to access them with "at" operators
        LSExpression taskArray = model.array();
        for (int j = 0; j < nbJobs; ++j) {
            taskArray.addOperand(model.array(tasks[j].begin(), tasks[j].end()));
        }

        LSExpression 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) {
            LSExpression intensityLambda =
                model.createLambdaFunction([&](LSExpression 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
            LSExpression sequence = jobsOrder[m];
            model.constraint(model.eq(model.count(sequence), nbJobs));

            // Disjunctive resource constraints between the tasks on a machine
            LSExpression sequenceLambda = model.createLambdaFunction([&](LSExpression 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 solver
        localsolver.getParam().setTimeLimit(timeLimit);

        localsolver.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) {
            LSCollection 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 %LS_HOME%\bin\localsolvernet.dll .
csc JobshopIntensity.cs /reference:localsolvernet.dll
JobshopIntensity instances\i01_ft06.txt
using System;
using System.IO;
using localsolver;

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;

    // LocalSolver
    private LocalSolver localsolver;

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

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

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

    public JobshopIntensity()
    {
        localsolver = new LocalSolver();
    }

    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()
    {
        localsolver.Dispose();
    }

    public void Solve(int timeLimit)
    {
        // Declare the optimization model
        LSModel model = localsolver.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 LSExpression[nbJobs, nbMachines];
        for (int j = 0; j < nbJobs; ++j)
        {
            for (int m = 0; m < nbMachines; ++m)
            {
                tasks[j, m] = model.Interval(0, timeHorizon);
            }
        }

        // Create LocalSolver arrays to be able to access them with "at" operators
        LSExpression taskArray = model.Array(tasks);
        LSExpression 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)
        {
            LSExpression 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 LSExpression[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
            LSExpression sequence = jobsOrder[m];
            model.Constraint(model.Count(sequence) == nbJobs);

            // Disjunctive resource constraints between the tasks on a machine
            LSExpression 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 solver
        localsolver.GetParam().SetTimeLimit(timeLimit);

        localsolver.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)
            {
                LSCollection 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 %LS_HOME%\bin\localsolver.jar
java -cp %LS_HOME%\bin\localsolver.jar;. JobshopIntensity instances\i01_ft06.txt
Compilation / Execution (Linux)
javac JobshopIntensity.java -cp /opt/localsolver_13_0/bin/localsolver.jar
java -cp /opt/localsolver_13_0/bin/localsolver.jar:. JobshopIntensity instances/i01_ft06.txt
import java.util.*;
import java.io.*;
import localsolver.*;

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;

    // LocalSolver
    final LocalSolver localsolver;
    // Decision variables: time range of each task
    private LSExpression[][] tasks;
    // Decision variables: sequence of tasks on each machine
    private LSExpression[] jobsOrder;
    // Objective = minimize the makespan: end of the last task of the last job
    private LSExpression makespan;

    public JobshopIntensity(LocalSolver localsolver, String fileName) throws IOException {
        this.localsolver = localsolver;
    }

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

        // Create LocalSolver arrays to be able to access them with "at" operators
        LSExpression taskArray = model.array(tasks);
        LSExpression 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) {
            LSExpression mExpr = model.createConstant(m);
            LSExpression intensityLambda = model.lambdaFunction(
                    t -> model.at(intensityArray, mExpr, t));
            for (int j = 0; j < nbJobs; ++j) {
                LSExpression 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 LSExpression[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
            LSExpression sequence = jobsOrder[m];
            model.constraint(model.eq(model.count(sequence), nbJobs));

            // Disjunctive resource constraints between the tasks on a machine
            LSExpression mExpr = model.createConstant(m);
            LSExpression 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 solver
        localsolver.getParam().setTimeLimit(timeLimit);

        localsolver.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) {
                LSCollection 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 (LocalSolver localsolver = new LocalSolver()) {
            JobshopIntensity model = new JobshopIntensity(localsolver, 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);
        }
    }
}