Our model is working but infeasible !

hello everyone, 

This model that we develop is working but infeasbile. We can't finde the probleme where it is!
it is in the code or it is in the constraints or in the number that we used ?
thank you so much 

this is the code in python

import sys
import gurobipy as gp
from gurobipy import GRB 

# Parameters

N= 6                             #tüm görevlerin kümesi
n=[0, 2, 2, 4, 4, 4]             #i görevi için gerçekleşecek kasa besleme sayısı  
   
R= 16                              #maksimum rota sayısı 
Q= 3                              #trenin taşıyabileceği maksimum kasa sayısı
M = 2000000                        #büyük bir sayı

u = [0, 48, 150, 40, 40, 40]         #maksimum parça seviyesi
v = [0, 3, 5, 2, 2, 10]              #minimum parça seviyesi


c = [0, 74, 36, 126, 126, 1]         #parça çevrim süresi (sn)

e =[[0, 0, 0, 0],
    [3552, 7104],                   #H2 in 1. referansı
    [5400, 10800],                  #H2 in 2. referansı
    [5040, 10080, 15120, 20160],    #H4 in 1. referansı     #c*u
    [5040, 10080, 15120, 20160],    #H4 in 2. referansı
    [40, 80, 120, 160]]             #H4 nin 3. referansı
#k. i görevi için kasanın en erken teslim zamanı
           


d =[[0, 0, 0, 0],
    [12600, 25200],                   #H2 in 1. referansı
    [12600, 25200],                   #H2 in 2. referansı
    [6300, 12600, 18900, 25200],      #H4 in 1. referansı
    [6300, 12600, 18900, 25200],      #H4 in 2. referansı
    [6300, 12600, 18900, 25200]]      #H4 nin 3. referansı
#k. i görevi için kasanın en geç teslim zamanı

p= [0, 3330, 5220, 4788, 4788, 30] #i görevinin süresi

w= [240, 15, 15, 15, 35, 35]      #i görevinin kasa başına çalışma süresi 

t= [[0, 45, 45, 45, 45, 55],
    [45, 0,  0, 15, 15, 25],
    [45, 0,  0, 15, 15, 25],
    [45, 15, 15, 0,  0, 10],
    [45, 15, 15, 0,  0, 10],
    [55, 25, 15, 10, 10, 0]]

#trenin görev i konumundan görev j konumuna seyahat süresi

#Create a new model
model = gp.Model("csknz")

# Variables

x = model.addVars(N,4,R,N,4, vtype=GRB.BINARY, name="x")
y = model.addVars(N,4, name="y")
s = model.addVars(N,4, name="s")


for j in range(N):
    for l in range(1, n[j]):
        for r in range(R):
         for i in range(N):
            for k in range(1, n[i]):
                    x[j,l,r,i,k] = model.addVar(vtype=GRB.BINARY)

                  

for i in range(N):
    for k in range(1, n[i]):
        s[i,k] = model.addVar()
       
        
for j in range(N):
    for l in range(1, n[j]):
        y[j,l] = model.addVar()
       
        
obj = gp.quicksum(t[i][j] * x[j,l,r,i,k] 
                  for i in range(N) 
                  for k in range(1, n[i])
                  for j in range(N) 
                  for l in range(1, n[j]) 
                  for r in range(R))
model.setObjective(obj, GRB.MINIMIZE)

#addConstraints
# for i in range(N):
#     if i != 0:
#         model.addConstr(p[i] == (u[i]-v[i]) * c[i]) # cons.2

# for i in range(N):
#     if i != 0:
#         for k in range(1, n[i]):
#               model.addConstr(e[i][k]+1 == e[i][k] + p[i]) # cons.3          

# for i in range(N):
#       if i != 0:
#           for k in range(1, n[i]):
#               model.addConstr(d[i][k] == e[i][k] + (v[i]- 0) * c[i]) # cons.4

# model.addConstr(x[3,1,0,0,1]==1)
# model.addConstr(x[4,1,0,3,1]==1) 

for i in range(N):
      if i != 0:
          for k in range(1, n[i]):
              model.addConstr(e[i][k] <= s[i,k])     # cons.5
              model.addConstr(s[i,k] <= d[i][k])
            
model.addConstr(gp.quicksum(x[j,l,r,i,k] for j in range(N) if j != 0
                                        for l in range(1,n[j]))==1)  # cons.6




model.addConstr(gp.quicksum(x[j,l,r,i,k]for j in range(N) if j != 0 
                                          for l in range(1,n[j]) 
                                          for r in range(R))<=1)  # cons.7

model.addConstrs((gp.quicksum(x[i,k,r,i,k] for i in range(N)
                              for k in range(1,n[i])) == 0 
                              for r in range(R))) # cons.8


for i in range(N):
    if i != 0:
      for k in range(1,n[i]):
        model.addConstr((gp.quicksum(x[j,l,r,i,k] for j in range(N) 
                              for l in range(1,n[j]) 
                              for r in range(R)) == 1))  # cons.9


for j in range(N):
    if j != 0 :
      for l in range(1,n[j]):
        model.addConstr((gp.quicksum(x[j,l,r,i,k] for i in range(N)
                                for k in range(1,n[i])  
                                for r in range(R)) == 1)) # cons.10




model.addConstrs((gp.quicksum(x[j,l,r,i,k]for i in range(N) 
                                  for k in range(1,n[i]) 
                                  for j in range(N) if j != 0 
                                  for l in range(1,n[j]))<= Q for r in range(R))) #cons.11
for i in range(N):
    for j in range(N):
        for k in range(1,n[i]):
            for l in range(1,n[j]):
                model.addConstr(s[i,k]+(w[i]+t[i][j]*gp.quicksum(x[j,l,r,i,k]for r in range(R)))-
        M*(1-gp.quicksum(x[j,l,r,i,k] for r in range(R)))+(y[j,l]-y[i,k])*(t[i][0]+w[0]+t[0][j]-t[i][j]) <= s[j,l])
#  #cons.12


for i in range(N):
    if i != 0 :              
      for l in range(1,n[j]):
        model.addConstr((gp.quicksum(r*x[j,l,r,i,k] for i in range(N) 
                                for k in range(1,n[i]) 
                                for r in range(R))==y[j,l]))  #cons.13

for i in range(N):
      for j in range(N):
        for k in range(1,n[i]):
          for l in range(1,n[j]):
              model.addConstr(y[j,l] >= (y[i,k]- M*(1-gp.quicksum(x[j,l,r,i,k] for r in range(R))))) #cons.14
                    
model.optimize()
status=model.status