CLSP with rework option and minimum lot size constraint modeling

Hi,

I tried to implement  CLSP-Rmin model from A.Goerler and S.Voss science article to Jupiter with Gurobi optimizer. 

When I applied CLSP without rework part (proportion of detective items O=0 instead of 0.01) it worked fine and gave correct solution. In this case the model does not use any parameters with 'r' but when I applied full CLSP with rework it shows me, that the model is infeasible, although in the research paper the results are Z=13560,25. 

What could be the cause that the model doesn't work? Is there a mistake in my formulation of the model?

Thank you in advance for help :)

#CLSP-Rmin example
from gurobipy import *
import time

start = time.time()

m = Model("CLSP-Rmin")
#Parameters
d= [[258, 248, 193, 213, 238, 276, 329, 359, 375, 362],
    [341, 322, 251, 288, 330, 376, 443, 476, 498, 471],
    [510, 473, 390, 422, 482, 570, 667, 712, 753, 721],
    [508, 480, 390, 418, 478, 580, 653, 726, 744, 712]]
f=   [110, 120, 100, 120]
h= [1.5, 2, 2.5, 1.25]
hR= [1.5, 2, 2.5, 1.25]
tp= [1, 1, 1, 1]
tpR= [0.75, 0.75, 0.75, 0.75]
ts= [25, 20, 20, 15]
tsR= [25, 20, 20, 15]
b= [2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000]
k= [0, 0, 0, 0]
#O= [0,0,0,0,0,0,0,0,0,0]
O= [0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01] #proportion of defective items
M= 100000
s0 = [0,0,0,0]
sR0 = [0,0,0,0]
T=len(b)

#Range of products and periods
Pro= range(len(f))
Pe = range(len(b))

#Decision variables default settings are 0.0 lb for all continous variables
x = m.addVars(Pro, Pe, vtype=GRB.BINARY, name="x")
s = m.addVars(Pro, Pe, vtype=GRB.INTEGER, name="s")
y = m.addVars(Pro, Pe, vtype=GRB.INTEGER, name="y")
yS = m.addVars(Pro, Pe, vtype=GRB.INTEGER, name="yS")
yR = m.addVars(Pro, Pe, vtype=GRB.INTEGER, name="yR")
sR = m.addVars(Pro, Pe, vtype=GRB.INTEGER, name="sR")
xR = m.addVars(Pro, Pe, vtype=GRB.BINARY, name="xR")
R = m.addVars(Pro, Pe, vtype=GRB.INTEGER, name="R")
m.update()

#Objective
obj= quicksum(((f[j]*x[j,t]+h[j]*s[j,t])+(f[j]*xR[j,t]+hR[j]*sR[j,t])) for j in Pro for t in Pe)
m.setObjective(obj, GRB.MINIMIZE) 
m.update()

#Capacity constraints
m.addConstrs(
    (quicksum(((tp[j]*y[j,t]+ts[j]*x[j,t])+(tpR[j]*yR[j,t]+tsR[j]*xR[j,t])) for j in Pro) <= b[t]
        for t in Pe), "Capacity-Ct")

#!production amount at least as size of lot size
for j in Pro:
    for t in Pe:
        m.addConstr(y[j,t]>=k[j]*x[j,t])
        m.addConstr(yS[j,t] ==(y[j,t]*(1-O[t])))
        m.addConstr(R[j,t]==(y[j,t]*O[t]))
    
#Inventory balance constraints
for j in Pro:
    for t in Pe:
        if t==0:
            m.addConstr(s[j,t]==0)#s0[j]+yS[j,t]+yR[j,t]-d[j][t])
        else:
            m.addConstr(s[j,t]==s[j,t-1]+yS[j,t]+yR[j,t]-d[j][t])
m.update()
for j in Pro:
    for t in Pe:
        if t==0:
            m.addConstr(sR[j,t]==0)#sR0[j]+R[j,t]-yR[j,t])
        else:
            m.addConstr(sR[j,t]==sR[j,t-1]+R[j,t]-yR[j,t])

m.update()
#s[j,T] the last period of product j equals 0
m.addConstrs((s[j,T-1] == 0) for j in Pro)    #T is defined as T=len(b) to count how many periods we have =10                 
m.addConstrs((sR[j,T-1] == 0) for j in Pro)
           
#Big M+ Big M for rework
m.addConstrs((y[j,t]<= sum(d[j][t:])*x[j,t] for j in Pro for t in Pe), "BigM-Ct")
m.addConstrs((yR[j,t]<= sum(d[j][t:])*xR[j,t] for j in Pro for t in Pe), "BigM-Ct for Rework")

logging.basicConfig(filename='info.log', level=logging.INFO)

m.optimize()
#m.computeIIS()

#m.write('model.ilp')
m.printAttr('x','y*')
m.printAttr('x','s*')
m.update()

#holding cost, cost for rework and setup cost chronologically
pp=0
for j in Pro:
    for t in Pe:
        pp= pp + (h[j]*s[j,t])
ppp=0
for j in Pro:
    for t in Pe:
        ppp= ppp + (hR[j]*sR[j,t])
pppp=0
for j in Pro:
    for t in Pe:
        pppp=pppp+(f[j]*x[j,t]+f[j]*xR[j,t])
        
print(pp.getValue())
print(ppp.getValue())
print(pppp.getValue())
#print(ppp)

and output like this: 

Explored 0 nodes (0 simplex iterations) in 0.01 seconds (0.00 work units)
Thread count was 1 (of 8 available processors)

Solution count 0

Model is infeasible
Best objective -, best bound -, gap -