Infeasible solution for Vehicle routing problem with outsourcing

Question

I have this Vehicle Routing Problem with the option to outsource part of your deliveries. Now, this code is based on a Mixed-Integer Linear Formulation made by a group of researchers and I cannot seem to find a mistake in my code that would cause the infeasibility. Anyone has some suggestions what could be the cause?

import pyomo.environ as pe
import pyomo.opt as po
import numpy as np
model = pe.ConcreteModel()

# Coordinates
xc = {0:40,1:25,2:22,3:22,4:20,5:20,6:18,7:15,8:15,9:10,10:10}
yc = {0:50,1:85,2:75,3:85,4:80,5:85,6:75,7:75,8:80,9:35,10:40}

# Amount of nodes
Customers = 10

# Customers -> NEW
N = {i for i in range (1, Customers+1)}

# Amount of nodes (SET) -> NEW
V = {0} | N  

# Set of arcs (SET)
A = {(i,j) for i in V for j in V if i!=j} 

# Demand
d = {0:0, 1:20, 2:30, 3:40, 4:10, 5:20, 6:20, 7:20, 8:10, 9:20, 10:30}
     
# Earliest service time 
e = {0:0, 1:673, 2:152, 3:471, 4:644, 5:73, 6:388, 7:300, 8:367, 9:371, 10:519}

# Latest service time
l = {0:960, 1:793, 2:272, 3:591, 4:764, 5:193, 6:508, 7:420, 8:487, 9:491, 10:639}

# Service time
s = {0:0, 1:10, 2:10, 3:10, 4:10, 5:10, 6:10, 7:10, 8:10, 9:10, 10:10}

# Travel time (from i to j) -> == eucledian distance between nodes based on coordinates
t = {(i,j): np.hypot(xc[i]-xc[j],yc[i]-yc[j]) for i,j in A} 

# Vehicle types (SET)
K = {1,2,3,4,5,6}

# Vehicle capacity
Q = {1:100, 2:200, 3:300, 4:400, 5:500, 6:1000}

# Fixed vehicle cost
f = {1:150, 2:350, 3:550, 4:800, 5:1100, 6:2500}

# Fleet size
m = {1:3 ,2:3 ,3:3 ,4:3 ,5:3 ,6:3}

# Travel cost for vehicle K (from i to j)
c = {(k,i,j): np.hypot(xc[i]-xc[j],yc[i]-yc[j]) for k in K for i,j in A} 

# Big M
M = 999999

# Cost function segments (SET)
S = {1,2,3,4}

# Lower bounds cost structure
L = {1:0, 2:200, 3:400, 4:600}

# Upper bounds cost structure
U = {1:200, 2:400, 3:600, 4:10000}

# Variable outsourcing costs
Ms = {1:3.50, 2:3.00, 3:2.50, 4:2.00}

# Fixed outsourcing costs
C = {1:100, 2:100, 3:100, 4:100}
##Sets

model.N = pe.Set(initialize = V, ordered = False) #Customers + Depot
model.N0= pe.Set(initialize = N, ordered = False) #Customers only
model.A = pe.Set(initialize = A, ordered = False) #Set of arcs
model.K = pe.Set(initialize = K, ordered = False) #Set of vehicles
model.S = pe.Set(initialize = S, ordered = False) #Set of discount segments
##Variables

# if a type−k vehicle traverses arc (i, j)
model.x = pe.Var(model.K, model.A, domain = pe.Binary) #mention model.A once because it is a combination of i and j

# vehicle’s load on arc (i, j)
model.y = pe.Var(model.A, domain = pe.NonNegativeReals) 

# vehicle’s arrival time at customer i
model.a = pe.Var(model.N, domain = pe.NonNegativeReals)

# if customer i is assigned to the common carrier
model.δ = pe.Var(model.N, domain = pe.Binary)

# quantity within the range [Ls, Us] of segment s to be subcontracted.
model.q = pe.Var(model.S, domain = pe.NonNegativeReals)

# if segment s is selected (different levels of fixed costs?)
model.z = pe.Var(model.S, domain = pe.Binary)
##Objective

ObjExpr1 = sum(f[k]*model.x[k,i,j] for k in model.K for (i,j) in model.A if i==0) #Fixed vehicle cost

ObjExpr2 = sum(c[k,i,j]*model.x[k,i,j] for k in model.K for (i,j) in model.A) #Travelling cost

ObjExpr3 = sum(Ms[s] * model.q[s] + C[s] * model.z[s] for s in model.S) #Outsourcing costs (note the Ms instead of M)

ObjExprfinal = ObjExpr1 + ObjExpr2 + ObjExpr3  
model.obj = pe.Objective(expr = ObjExprfinal, sense = pe.minimize)
##Constraints

# Visit a customer only once (CAUSE FOR INFEASIBILITY)
model.VisitOnce = pe.ConstraintList()
for j in model.N0:
    expression = sum(model.x[k,i,j] + model.δ[j] for k in model.K for (i,j) in model.A) == 1
    model.VisitOnce.add(expression)

# Flow conservation
model.FlowCons = pe.ConstraintList()
for k in model.K:
    for i in model.N:
        expression = sum(model.x[k,i,j] for (i,j) in model.A) - sum(model.x[k,j,i] for (j,i) in model.A) == 0
        model.FlowCons.add(expression)

# Vehicle availability
model.VehAvai = pe.ConstraintList()
for k in model.K:
    expression = sum(model.x[k,i,j] for (i,j) in model.A if i==0) <= m[k]
    model.VehAvai.add(expression)

# Demand satisfaction (CAUSE FOR INFEASIBILITY)
model.DemandSatisf = pe.ConstraintList()
for i in model.N0:
    expression = sum(model.y[i,j] for (i,j) in model.A) - sum(model.y[j,i] for (j,i) in model.A) == d[i]*(1 - model.δ[i])
    model.DemandSatisf.add(expression)

# Vehicle capacity
model.VehCap = pe.ConstraintList()
for k in model.K:
    for (i,j) in model.A:
        expression = model.y[i,j] <= sum((Q[k]-d[i])*model.x[k,i,j] for k in model.K)
        model.VehCap.add(expression)

# Planning must follow travel times
model.Planning = pe.ConstraintList()
for (i,j) in model.A:
    if i !=0:
        expression = model.a[j]-model.a[i] >= s[i]+t[i,j] - (M*(1-sum(model.x[k,i,j] for k in model.K)))
        model.Planning.add(expression)

# First customer may not be visited before...
model.FirstCust = pe.ConstraintList()
for j in model.N:
    if j !=0:
        expression = model.a[j] >= t[0,j]*sum(model.x[k,0,j] for k in model.K)
        model.FirstCust.add(expression)

# Cost related
model.cost9 = pe.ConstraintList()
expression = sum(d[i]*model.δ[i] for i in model.N0) <= sum(model.q[s] for s in model.S)
model.cost9.add(expression)

# Cost related  
model.cost10 = pe.ConstraintList()
for s in model.S:
    expression = (L[s]*model.z[s]) <= (model.q[s])
    model.cost10.add(expression)

# Cost related  
model.cost101 = pe.ConstraintList()
for s in model.S:
    expression = (model.q[s]) <= (U[s]*model.z[s])
    model.cost101.add(expression)    
    
# Cost related
model.cost11 = pe.ConstraintList()
expression = sum(model.z[s] for s in model.S) == 1
model.cost11.add(expression)

# Arrival must be in time window
model.TimeWind = pe.ConstraintList()
for i in model.N:
    expression = e[i] <= model.a[i] 
    model.TimeWind.add(expression)

# Arrival must be in time window
model.TimeWind1 = pe.ConstraintList()
for i in model.N:
    expression = model.a[i] <= l[i]
    model.TimeWind1.add(expression)
#model.print()
solver = po.SolverFactory('gurobi')
result = solver.solve(model, tee = True)
Set parameter Username
Academic license - for non-commercial use only - expires 2022-07-09
Read LP format model from file C:\Users\Danie\AppData\Local\Temp\tmpvfas1ihl.pyomo.lp
Reading time = 0.01 seconds
x801: 795 rows, 801 columns, 11526 nonzeros
Gurobi Optimizer version 9.5.0 build v9.5.0rc5 (win64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 795 rows, 801 columns and 11526 nonzeros
Model fingerprint: 0x66a8aa5e
Variable types: 126 continuous, 675 integer (675 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+04]
  Objective range  [2e+00, 3e+03]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+03]
Presolve removed 729 rows and 111 columns
Presolve time: 0.00s

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

Solution count 0

Model is infeasible
Best objective -, best bound -, gap -
WARNING: Loading a SolverResults object with a warning status into
    model.name="unknown";
      - termination condition: infeasible
      - message from solver: Model was proven to be infeasible.
print(result.solver.status)
print(result.solver.termination_condition)
warning
infeasible