PuLP Python on 3D bin packing problem / container loading problem

Question

I am trying to make a program with PuLP on the mathematical model in this journal: An analytical model for the container loading problem https://doi.org/10.1016/0377-2217(94)00002-T https://www.sciencedirect.com/science/article/pii/037722179400002T

Here's my code:

    #objective function
prob += lpSum(L[j]*W[j]*H[j] for j in range(m)) - lpSum(p[i]*q[i]*r[i] for i in range(N))
    #constraints

for i in range(N):
    for k in range(N):
      if i<k:
        prob.addConstraint(x[i]+p[i]*l_x[i]+q[i]*(l_z[i]-w_y[i]+h_z[i])+r[i]*(1-l_x[i]-l_z[i]+w_y[i]-h_z[i])<=x[k]+(1-a[(i,k)])*M)
      else:
        pass

for i in range(N):
    for k in range(N):
      if i<k:
        prob.addConstraint(x[k]+p[k]*l_x[k]+q[k]*(l_z[k]-w_y[k]+h_z[k])+r[k]*(1-l_x[k]-l_z[k]+w_y[k]-h_z[k])<=x[i]+(1-b[(i,k)])*M)
      else:
        pass

for i in range(N):
    for k in range(N):
      if i<k:
        prob.addConstraint(y[i]+q[i]*w_y[i]+p[i]*(1-l_x[i]-l_z[i])+r[i]*(l_x[i]+l_z[i]-w_y[i])<=y[k]+(1-c[(i,k)])*M)
      else:
        pass

for i in range(N):
    for k in range(N):
      if i<k:
        prob.addConstraint(y[k]+q[k]*w_y[k]+p[k]*(1-l_x[k]-l_z[k])+r[k]*(l_x[k]+l_z[k]-w_y[k])<=y[i]+(1-d[(i,k)])*M)
      else:
        pass

for i in range(N):
    for k in range(N):
      if i<k:
        prob.addConstraint(z[i]+r[i]*h_z[i]+q[i]*(1-l_z[i]-h_z[i])+p[i]*l_z[i]<=z[k]+(1-e[(i,k)])*M)
      else:
        pass

for i in range(N):
    for k in range(N):
      if i<k:
        prob.addConstraint(z[k]+r[k]*h_z[k]+q[k]*(1-l_z[k]-h_z[k])+p[k]*l_z[k]<=z[i]+(1-f[(i,k)])*M)
      else:
        pass

for i in range(N):
    for j in range(m):
        prob.addConstraint(x[i]+p[i]*l_x[i]+q[i]*(l_z[i]-w_y[i]+h_z[i])+r[i]*(1-l_x[i]-l_z[i]+w_y[i]-h_z[i])<=L[j]+(1-s[(i,j)])*M)

for i in range(N):
    for j in range(m):
        prob.addConstraint(y[i]+q[i]*w_y[i]+p[i]*(1-l_x[i]-l_z[i])+r[i]*(l_x[i]+l_z[i]-w_y[i])<=W[j]+(1-s[(i,j)])*M)

for i in range(N):
    for j in range(m):
        prob.addConstraint(z[i]+r[i]*h_z[i]+q[i]*(1-l_z[i]-h_z[i])+p[i]*l_z[i]<=H[j]+(1-s[(i,j)])*M)

for i in range(N):
  prob.addConstraint(l_x[i]+l_y[i]+l_z[i]==1)

for i in range(N):
  prob.addConstraint(w_x[i]+w_y[i]+w_z[i]==1)

for i in range(N):
  prob.addConstraint(h_x[i]+h_y[i]+h_z[i]==1)

for i in range(N):
  prob.addConstraint(l_x[i]+w_x[i]+h_x[i]==1)

for i in range(N):
  prob.addConstraint(l_y[i]+w_y[i]+h_y[i]==1)

for i in range(N):
  prob.addConstraint(l_z[i]+w_z[i]+h_z[i]==1)

for i in range(N):
  for j in range(m):
    for k in range(N):
      if i<k:
        prob.addConstraint(a[(i,k)]+b[(i,k)]+c[(i,k)]+d[(i,k)]+e[(i,k)]+f[(i,k)]>=s[(i,j)]+s[(k,j)]-1)
      else:
        pass

for i in range(N):
  prob.addConstraint(lpSum(s[(i,j)] for j in range(m))==1)

for j in range(m):
  prob.addConstraint(lpSum(s[(i,j)] for i in range(N))<=M*n[j])

When i tried to solve it by

prob.solve()
LpStatus[prob.status]

The output says that my problem is infeasible, and also the vast majority of the varValue is zero. can you help me fix my code? i think there's a mistake in declaring the constraints and obj. function.

thanks!