from gurobipy import *
import numpy as np
import csv

# Parameters needed are:
(1) the total number of jobs (n). Here I denote it by "NumofJobs"
(2) the total number of machines (m). Here I denote it by "NumofMachines"
(3) the processing times.  Here I use a numpy matrix: "Tasktime[r, j]" : p_{r,j}


# Building the model:
m = Model ("Wilson")


# Generating variables:
z = {}
for j in range(NumofJobs):
    for i in range(NumofJobs):
        z[j, i] = m.addVar(vtype=GRB.BINARY)

s = {}
for r in range(NumofMachines):
    for j in range(NumofJobs):
        s[r, j] = m.addVar(vtype=GRB.CONTINUOUS)

m.update()


# Generating constraints:
for i in range(NumofJobs):
    m.addConstr(quicksum(z[j, i] for j in range(NumofJobs)) == 1)

for j in range(NumofJobs):
    m.addConstr(quicksum(z[j, i] for i in range(NumofJobs)) == 1)

m.addConstr(s[0, 0] == 0)

for i in range(NumofJobs - 1):
    m.addConstr(s[0, i] + quicksum((Tasktime[0, j]) * z[j, i] for j in range(NumofJobs)) == s[0, i + 1])

for r in range(NumofMachines - 1):
    m.addConstr(s[r, 0] + quicksum((Tasktime[r, j]) * z[j, 0] for j in range(NumofJobs)) == s[r + 1, 0])

for r in range(NumofMachines - 1):
    for i in range(1, NumofJobs):
        m.addConstr(s[r, i] + quicksum(Tasktime[r, j] * z[j, i] for j in range(NumofJobs)) <= s[r + 1, i])

for r in range(1, NumofMachines):
    for i in range(NumofJobs - 1):
        m.addConstr(s[r, i] + quicksum(Tasktime[r, j] * z[j, i] for j in range(NumofJobs)) <= s[r, i + 1])


# Setting the objective function:
m.setObjective(s[NumofMachines - 1, NumofJobs - 1] + quicksum(Tasktime[NumofMachines - 1, j] * z[j, NumofJobs - 1] for j in range(NumofJobs)), GRB.MINIMIZE)


m.optimize()