<jupyter_start><jupyter_text>Tutorial 4: Scheduling with CP: Some Basics  Today we discover a new family of problems, called scheduling. Scheduling problems are widely present in real life applications such as timetabelling, transportation, project management, and manufacturing. We consider a particular problem called the job shop scheduling problem.  In this problem, we are given a set of $n$ jobs: $J_1, J_2, \ldots,  J_n$ and a set of $m$ machines $M_1, M_2, \ldots,  M_m$. Each job $J_i$ is defined by a set of $m$ (non-preemptive) tasks $T_{i,1} \ldots T_{i,m}$. Every task $T_{i,k}$ is associated with a duration $p_{i,k}$ and is supposed to be scheduled on machine $M_k$. The problem has two sets of constraints:  - Precedence constraints: Each job is associated with an order of tasks to respect when scheduling the different tasks. - Disjunctive constraints: Each machine can process only one task at a given timeThe standard optimisation version of this problem asks to minimize the makespan, i.e., the total scheduling time. Constraint programming has been widely and successfully used to solve scheduling problems. Many global constraints have been proposed. CP solvers often offer a dedicated library for scheduling. Please have a look at the diffrent features proposed in docplex here http://ibmdecisionoptimization.github.io/docplex-doc/cp/docplex.cp.modeler.py.html?highlight=schedulingscheduling-functions  In this tutorial, we will be using (only):  - Interval variables for the different variables of the problem:  http://ibmdecisionoptimization.github.io/docplex-doc/cp/docplex.cp.expression.py.html?highlight=interval_vardocplex.cp.expression.interval_var  - end_before_start constraints to model precedence constraints : http://ibmdecisionoptimization.github.io/docplex-doc/cp/docplex.cp.modeler.py.html?highlight=end_before_startdocplex.cp.modeler.end_before_start- no_overlap global constraint : http://ibmdecisionoptimization.github.io/docplex-doc/cp/docplex.cp.modeler.py.html?highlight=no_overlapdocplex.cp.modeler.no_overlap The format for a job shop scheduling instance respects the following syntax:  - The first line containts only $n$ $m$ in this order  ($n$ is the number of jobs and $m$ is the number of machines) - Then $n$ lines are given. Each line $i$ is associated to the job $J_i$ and contains exactly $2 m$ integers $x^i_1$ $d^i_1$ $x^i_2$ $d^i_2$ $\ldots$ $x^i_n$ $d^i_n$. Each $x^i_k$ is the $k^{th}$ machine required by the $k^{th}$ task of the job $J_i$, and $d^i_k$ represents its duration.  Consider for example the data file instance.data :      6  6     2   1  0   3  1   6  3   7  5   3  4   61   8  2   5  4  10  5  10  0  10  3   42   5  3   4  5   8  0   9  1   1  4   7  1   5  0   5  2   5  3   3  4   8  5   9  2   9  1   3  4   5  5   4  0   3  3   1  1   3  3   3  5   9  0  10  4   4  2   1 This instance has $6$ jobs and $6$ machines. The first job requires the execution of task $T_{1,2}$ (on machine $2$) of duration $1$, then $T_{1,0}$ (on machine $0$) of duration $3$, etc.  Write a simple python code to parse the data file instance.data<jupyter_code>import numpy as np
datapath = './example.data'

def read_data(datapath):
    jobs = []

    with open(datapath,'r') as fp:
        for line in fp:
            tasks = []
            if len(line.split()) == 2: 
                NB_MACHINES = int(line.split()[1])
            else :
                tasks = list(map(lambda x : int(x),line.split()))
                tasks = np.array(tasks).reshape(NB_MACHINES,2)
                jobs.append(tasks)    
    return np.array(jobs)

jobs = read_data(datapath)
NB_JOBS = jobs.shape[0]
NB_MACHINES = jobs.shape[1]
print("Nb of jobs : ", NB_JOBS)
print("Nb of machines : ", NB_MACHINES)
print("First job requires the execution of :")
for i in range(NB_MACHINES):
    print("Task %d on machine %d of duration %d" % (i,jobs[0][i][0],jobs[0][i][1]))<jupyter_output>First job requires the execution of :
Task 0 on machine 2 of duration 1
Task 1 on machine 0 of duration 3
Task 2 on machine 1 of duration 6
Task 3 on machine 3 of duration 7
Task 4 on machine 5 of duration 3
Task 5 on machine 4 of duration 6<jupyter_text>Create a matrix machine that satisfies: machine[i][k] is the machine of the $k^{th}$ task of job $i$ (one line)<jupyter_code>machine = np.ones((NB_JOBS, NB_MACHINES))                       # initialize with ones
for id_job in range(NB_JOBS):
    for id_task in range(NB_MACHINES):
        machine[id_job][id_task] = int(jobs[id_job][id_task][0])
print(machine)<jupyter_output>[[2. 0. 1. 3. 5. 4.]
 [1. 2. 4. 5. 0. 3.]
 [2. 3. 5. 0. 1. 4.]
 [1. 0. 2. 3. 4. 5.]
 [2. 1. 4. 5. 0. 3.]
 [1. 3. 5. 0. 4. 2.]]<jupyter_text>Create a matrix duration that satisfies: duration[i][k] is the duration of the $k^{th}$ task of job $i$ (one line)<jupyter_code>duration = np.zeros((NB_JOBS, NB_MACHINES))                       # initialize with zeros
for id_job in range(NB_JOBS):
    for id_task in range(NB_MACHINES):
        duration[id_job,id_task] = int(jobs[id_job,id_task,1])
print(duration)<jupyter_output>[[ 1.  3.  6.  7.  3.  6.]
 [ 8.  5. 10. 10. 10.  4.]
 [ 5.  4.  8.  9.  1.  7.]
 [ 5.  5.  5.  3.  8.  9.]
 [ 9.  3.  5.  4.  3.  1.]
 [ 3.  3.  9. 10.  4.  1.]]<jupyter_text>Compute a naive upper bound for the makespan. Note: this upper bound will be used as an upper bound for every interval variable we create.<jupyter_code>naive_ub = np.sum(np.sum(duration))
print("A naive upper bound = %d, where we respectively finish each task of each jobs" % naive_ub)<jupyter_output>A naive upper bound = 197, where we respectively finish each task of each jobs<jupyter_text>Create a CpoModel() and the different interval variables you need to solve this problem (don't forget to use the upper bound you computed earlier)<jupyter_code>from docplex.cp.model import CpoModel
from docplex.cp.modeler import all_diff
from config import setup
setup()

model = CpoModel(name="Job scheduling problem")
job_operations = [[model.interval_var(end=(0,int(naive_ub)), size=int(duration[j][m]), name="O{}-{}".format(j, m)) 
                   for m in range(NB_MACHINES)] for j in range(NB_JOBS)]<jupyter_output>/Users/msiala/.local/lib/python3.7/site-packages/docplex/cp/__init__.py:30: RuntimeWarning: DOcplex.CP is supported by Python versions 2.7.9+, 3.4.x, 3.5.x and 3.6.x, not '3.7.6.final.0'. Use it at your own risk.
  warnings.warn(msg, RuntimeWarning)<jupyter_text>Post the precedence constraints using the end_before_start constraints<jupyter_code>for id_job in range(NB_JOBS):
    for id_task in range(NB_MACHINES-1):
        model.add(model.end_before_start(job_operations[id_job][id_task],job_operations[id_job][id_task+1]))<jupyter_output><empty_output><jupyter_text>Post the disjunctive constraints using the no-overlap constraints<jupyter_code>machine_op = [[] for i in range(NB_JOBS)]
for id_job in range(NB_JOBS):
    for id_task in range(NB_MACHINES):
        print(id_job, id_task)
        machine_op[int(machine[id_job][id_task])].append(job_operations[id_job][id_task]) 
for id_machine in range(NB_MACHINES):
    model.add(model.no_overlap(machine_op[id_machine]))<jupyter_output>0 0
0 1
0 2
0 3
0 4
0 5
1 0
1 1
1 2
1 3
1 4
1 5
2 0
2 1
2 2
2 3
2 4
2 5
3 0
3 1
3 2
3 3
3 4
3 5
4 0
4 1
4 2
4 3
4 4
4 5
5 0
5 1
5 2
5 3
5 4
5 5<jupyter_text>Create a makespan interval variable and link it to some variables using the max global constrainthttp://ibmdecisionoptimization.github.io/docplex-doc/cp/docplex.cp.modeler.py.html?highlight=maxdocplex.cp.modeler.max<jupyter_code>max_global_constraint = model.max([model.end_of(job_operations[i][NB_MACHINES - 1]) for i in range(NB_JOBS)])<jupyter_output><empty_output><jupyter_text>Add now the makespan as an objective function<jupyter_code>model.add(model.minimize(max_global_constraint))<jupyter_output><empty_output><jupyter_text>Solve this instance. What is the value of the makespan you found. You can print the solution in a format that is easy to see visually.<jupyter_code>model.print_information()
sol = model.solve() #, SearchType='Restart',DefaultInferenceLevel = "Medium")
sol.print_solution()<jupyter_output>-------------------------------------------------------------------------------
Model constraints: 36, variables: integer: 0, interval: 36, sequence: 6
Solve status: Optimal
Search status: SearchCompleted, stop cause: SearchHasNotBeenStopped
Solve time: 0.01 sec
-------------------------------------------------------------------------------
Objective values: (55,)
          bounds: (55,)
          gaps: (0,)
O0-0: (start=5, end=6, size=1, length=1)
O0-1: (start=6, end=9, size=3, length=3)
O0-2: (start=16, end=22, size=6, length=6)
O0-3: (start=30, end=37, size=7, length=7)
O0-4: (start=42, end=45, size=3, length=3)
O0-5: (start=49, end=55, size=6, length=6)
O1-0: (start=0, end=8, size=8, length=8)
O1-1: (start=8, end=13, size=5, length=5)
O1-2: (start=13, end=23, size=10, length=10)
O1-3: (start=28, end=38, size=10, length=10)
O1-4: (start=38, end=48, size=10, length=10)
O1-5: (start=48, end=52, size=4, length=4)
O2-0: (start=0, end=5, size=5, length=5)
O2-1: (start=5, end=9, size=4, l[...]<jupyter_text>Factorise your code so that it takes as input the data file and it solves the problem. Try to use few other instances https://github.com/tamy0612/JSPLIB/tree/master/instances<jupyter_code>def solve_jobs_scheduling (datapath):
    # ========================
    # 
    # ========================
    jobs = read_data(datapath)
    NB_JOBS = jobs.shape[0]
    NB_MACHINES = jobs.shape[1]
    print("Nb of jobs : ", NB_JOBS)
    print("Nb of machines : ", NB_MACHINES)
    
    # ========================
    # Create a matrix of machine for each task
    # ========================
    machine = np.ones((NB_JOBS, NB_MACHINES))                       # initialize with ones
    for id_job in range(NB_JOBS):
        for id_task in range(NB_MACHINES):
            machine[id_job][id_task] = int(jobs[id_job][id_task][0])
    print(machine)
    # ========================
    # Create a matrix of duration for each task
    # ========================
    duration = np.zeros((NB_JOBS, NB_MACHINES))                       # initialize with zeros
    for id_job in range(NB_JOBS):
        for id_task in range(NB_MACHINES):
            duration[id_job,id_task] = int(jobs[id_job,id_task,1])
    
    
    # ========================
    # Compute the naive upper bound
    # ========================
    naive_ub = np.sum(np.sum(duration))
    print("A naive upper bound = %d, where we respectively finish each task of each jobs" % naive_ub)
    
    # ========================
    # Create CpoModel & Initialize interval variables for each task
    # ========================
    model = CpoModel(name="Job scheduling problem")
    job_operations = [[model.interval_var(end=(0,int(naive_ub)), size=int(duration[j][m]), name="O{}-{}".format(j, m)) 
                   for m in range(NB_MACHINES)] for j in range(NB_JOBS)]
    
    
    # ========================
    # Add Precedence constraints
    # ========================
    for id_job in range(NB_JOBS):
        for id_task in range(NB_MACHINES-1):
            model.add(model.end_before_start(job_operations[id_job][id_task],job_operations[id_job][id_task+1]))
            
    # ========================
    # Add Disjunctive constraints
    # ========================
    machine_op = [[] for i in range(NB_JOBS)]
    for id_job in range(NB_JOBS):
        for id_task in range(NB_MACHINES):
            machine_op[int(machine[id_job][id_task])].append(job_operations[id_job][id_task]) 
    for id_machine in range(NB_MACHINES):
        model.add(model.no_overlap(machine_op[id_machine]))
    
    # ========================
    # Add objective function
    # ========================
    max_global_constraint = model.max([model.end_of(job_operations[i][NB_MACHINES - 1]) for i in range(NB_JOBS)])
    model.add(model.minimize(max_global_constraint))
    
    # ========================
    # Print information of created model 
    # ========================
    model.print_information()
    
    # ========================
    # Solve the problem
    # ========================
    sol = model.solve() #, SearchType='Restart',DefaultInferenceLevel = "Medium")
    sol.print_solution()
    
    return sol
msol = solve_jobs_scheduling("example.data")<jupyter_output>Nb of jobs :  6
Nb of machines :  6
[[2. 0. 1. 3. 5. 4.]
 [1. 2. 4. 5. 0. 3.]
 [2. 3. 5. 0. 1. 4.]
 [1. 0. 2. 3. 4. 5.]
 [2. 1. 4. 5. 0. 3.]
 [1. 3. 5. 0. 4. 2.]]
A naive upper bound = 197, where we respectively finish each task of each jobs
Model: Job scheduling problem
 - expressions: 37
 - modeling time: 0.03 sec
-------------------------------------------------------------------------------
Model constraints: 36, variables: integer: 0, interval: 36, sequence: 6
Solve status: Optimal
Search status: SearchCompleted, stop cause: SearchHasNotBeenStopped
Solve time: 0.01 sec
-------------------------------------------------------------------------------
Objective values: (55,)
          bounds: (55,)
          gaps: (0,)
O0-0: (start=5, end=6, size=1, length=1)
O0-1: (start=6, end=9, size=3, length=3)
O0-2: (start=16, end=22, size=6, length=6)
O0-3: (start=30, end=37, size=7, length=7)
O0-4: (start=42, end=45, size=3, length=3)
O0-5: (start=49, end=55, size=6, length=6)
O1-0: (start=0, [...]<jupyter_text>At this stage you are completly free to play. Try different instances, different configurations of the solver, different branching strategies, restarts, randomisation, etc. You may want to present your results as a cactus 🌵 plot : in the x-axis we have the runtime, on the y-axis, we have the number of instances solved. Also, some instances are still open in the literature. Have a look here for an up to date list of bounds http://optimizizer.com/TA.php<jupyter_code>#Ceci permet d'afficher les solutions. Je n'arrive pas à débugger.. 

import docplex.cp.utils_visu as visu
print (machine)

if msol and visu.is_visu_enabled():
    visu.timeline("Solution for job-shop " + datapath)
    visu.panel("Jobs")
    for i in range(NB_JOBS):
        visu.sequence(name='J' + str(i),
                      intervals=[(msol.get_var_solution(job_operations[i][j]), machine[i][j], 'M' + str(machine[i][j])) for j in
                                 range(NB_MACHINES)])
    visu.panel("Machines")
    for k in range(NB_MACHINES):
        visu.sequence(name='M' + str(k),
                      intervals=[(msol.get_var_solution(machine_operations[k][i]), k, 'J' + str(i)) for i in range(NB_JOBS)])
    visu.show()<jupyter_output>[[2. 0. 1. 3. 5. 4.]
 [1. 2. 4. 5. 0. 3.]
 [2. 3. 5. 0. 1. 4.]
 [1. 0. 2. 3. 4. 5.]
 [2. 1. 4. 5. 0. 3.]
 [1. 3. 5. 0. 4. 2.]]<jupyter_text>What did you learn loday?<jupyter_code>We learnt how to resolve scheduling problem with CSP, by considering the two following global constraints : the end_before_start and no_overlap.<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16734612043039498, 'fraction_numerical': 0.03842812270266658, 'mean_word_length': 3.272986864648772, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 44, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '166009', 'license': 'CC-BY-3.0', 'n_tokens_mistral': 5545, 'n_tokens_neox': 5087, 'n_words': 1598, 'path': 'teaching/2021-2022/5-cp/Tutorial4-correction.ipynb', 'repo_name': 'siala/siala.github.io', 'size': '14963'}