The model is infeasible

Hello, 

I am trying to solve Unrelated Parallel Machine Scheduling with Release Time, 
and I am quoting the mathematical model from a published paper. Because it was published, I believe it is the right model, but it is returning infeasible.
Can you please help me to solve this issue? I will give the code in here. 

import numpy as np
import gurobipy as gp
from gurobipy import GRB

# number of jobs
n = 5

# number of machines
m = 3

# maximum time horizon
Tmax = 10

# processing time of job i on machine k
P = np.random.randint(1, 5, size=(n, m))

# release time of job i
r = np.random.randint(0, 3, size=n)

# big M
M = 1000

# Create a new Gurobi model
model = gp.Model("Parallel Machine Scheduling")

# Create decision variables
X = model.addVars(range(n), range(m), vtype=GRB.BINARY, name="X")
Y = model.addVars(range(n), range(n), range(m), vtype=GRB.BINARY, name="Y")
Z = model.addVars(range(n), range(m), range(Tmax), vtype=GRB.BINARY, name="Z")
f = model.addVars(range(n), lb=0, name="f")
t = model.addVars(range(n), lb=0, name="t")
Cmax = model.addVar(lb=0, name="Cmax")

# Set objective function
model.setObjective(Cmax, GRB.MINIMIZE)

# Add constraints

# 1. Calculates maximum completion time of all jobs
for i in range(n):
    model.addConstr(f[i] - Cmax <= 0)

# 2. Ensure job assigned to exactly one machine
for i in range(n):
    model.addConstr(sum(X[i, k] for k in range(m)) == 1)

for k in range(m):
    model.addConstr(sum(X[i, k] for i in range(n)) == 1)

# 3. To choose one binary variable for the setup time for each job on one machine with considering the immediate predecessor job
for j in range(n):
    for k in range(m):
        model.addConstr(sum(Y[q, j, k] for q in range(n) if q != j) == X[j, k])

for i in range(n):
    for k in range(m):
        model.addConstr(sum(Y[i, j, k] for j in range(n) if i != j) == X[i, k])

for k in range(m):
    model.addConstr(sum(Y[0, j, k] for j in range(n)) == 1)

# 4. Calculations of completion jobs
for i in range(n):
    model.addConstr(f[i] == t[i] + sum(P[i, k] * X[i, k] for k in range(m)))

# 5. Assign dummy job = 0
model.addConstr(f[0] == 0)

# 6. Guarantee job cannot start before its predecessor job
for j in range(n):
    for q in range(n):
        model.addConstr(f[q] - t[j] + M * (sum(Y[q, j, k] for k in range(m)) - 1) <= 0)

# 7. Release time constraint
for i in range(n):
    model.addConstr(r[i] - t[i] <= 0)

# 8. Time index constraint
for i in range(n):
    model.addConstr(sum(Z[i, k, T] for k in range(m) for T in range(Tmax)) == f[i] - t[i])

for i in range(n):
    for k in range(m):
        model.addConstr(sum(Z[i, k, T] for T in range(Tmax)) <= M * X[i, k])

for i in range(n):
    for T in range(Tmax):
        model.addConstr(f[i] >= T * sum(Z[i, k, T] for k in range(m)))
        model.addConstr(t[i] <= T * sum(Z[i, k, T] for k in range(m)) + M * (1 - (sum(Z[i, k, T] for k in range(m)))))

# Optimize the model
model.optimize()