title: Infeasible solution for Vehicle routing problem with outsourcing
body: <p>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?</p>
  <pre><code>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 -&gt; NEW
  N = {i for i in range (1, Customers+1)}
  
  # Amount of nodes (SET) -&gt; 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) -&gt; == 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) &lt;= 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] &lt;= 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] &gt;= 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] &gt;= 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) &lt;= 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]) &lt;= (model.q[s])
      model.cost10.add(expression)
  
  # Cost related  
  model.cost101 = pe.ConstraintList()
  for s in model.S:
      expression = (model.q[s]) &lt;= (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] &lt;= 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] &lt;= 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=&quot;unknown&quot;;
        - termination condition: infeasible
        - message from solver: Model was proven to be infeasible.
  print(result.solver.status)
  print(result.solver.termination_condition)
  warning
  infeasible
  </code></pre>