Obtaining wrong solutions for VRPTW with SDLs

Dear Gurobi community,

currently I am trying to implement a VRP with (hard) time windows extended by the opportunity to deliver packages to Shared Delivery Locations (SDL) instead of delivering customers who are served by these SDLs.

My implementation does provide feasible solutions, but these solutions are obviously totally wrong. I tried a lot, but I couldn't find the reason for these wrong results.

This is the Python code of the implemented model:

import numpy as np
import gurobipy as gp
import random
import time

random.seed(10)
def VRPTW_SDL(N,N_D,I,F,V,e,Tmax,a,l,delta,gamma,t,c,C):
    start_time = time.time()
    
    mod = gp.Model()
    
    x = mod.addVars(N_D,N_D,vtype = 'B', name = 'x') #binary variable indicating whether j is visited directly after node i (1) or not (0)
    Y = mod.addVars(N,F,vtype = 'B', name = 'Y') #binary variable indicating whether i is delivered to SDL (Shared Delivery Location) f (1) or not (0)
    z = mod.addVars(F,vtype ='B', name = 'z') #binary variable indicating whether SDL f is visited (1) or not
    s = mod.addVars(N_D,vtype = 'C',name = 's') #non-negative variable indicating the visit time at node i
    
    
    mod.setObjective(gp.quicksum(c[i,j]*x[i,j] for i in N_D for j in N_D if j!=i)
                     + gp.quicksum(delta*Y[i,f] for i in I for f in F) 
                     + gp.quicksum(gamma*x[0,j] for j in N), gp.GRB.MINIMIZE)
    
    mod.addConstrs((gp.quicksum(x[i,j] for i in N_D if i!=j) 
                    + gp.quicksum(Y[i,f] for f in V[j-1]) 
                    == 1 for j in I),name = 'NB_1')
    
    mod.addConstrs((gp.quicksum(x[j,i] for j in N_D ) ==
                    gp.quicksum(x[i,j] for j in N_D ) for i in N_D), name = 'NB_2')
    
    mod.addConstrs((z[f] >= (1/len(I))*gp.quicksum(Y[i,f] for i in I) 
                                            for f in F), name = 'NB_3')
    
    mod.addConstrs((gp.quicksum(x[i,f] for i in N_D) == z[f] for f in F), name = 'NB_4')
    
    mod.addConstrs((s[i] + t[i,j] + e[i] - (2*Tmax)*(1-x[i,j]) <= s[j] for i in N_D for j in N if i!=j),name = 'NB_5')
    
    mod.addConstrs(((-Tmax*(gp.quicksum(Y[i,f] for f in F) + a[i])) <= s[i] for i in I), name= 'NB_6a')
    
    mod.addConstrs(((Tmax*(gp.quicksum(Y[i,f] for f in F) + l[i])) >= s[i] for i in I), name= 'NB_6b')
    
    mod.addConstrs((s[j] + e[j] + t[j,0] <= Tmax for j in N), name = 'NB_7')
    
    mod.addConstrs((gp.quicksum(Y[i,f] for i in I) <= B[f-C-1] for f in F), name = 'NB_8')
    
    mod.addConstrs((Y[i,f] == 0 for i in I for f in F if f not in V[i-1]), name = 'NB_9')
    
    mod.optimize()
    
    mod.write('VRPTW_SDL.mps')
    
    #If an optimal is found, print values of decision variables and objective function
    if mod.status == gp.GRB.OPTIMAL:
        print(mod.printAttr('X'))
        print(mod.ObjVal)
    #if model is infeasible, print those constraints that are not satisfied
    elif mod.status == gp.GRB.Status.INFEASIBLE:
        print('Optimization was stopped with status %d' % mod.status)
        mod.computeIIS()
        for c in mod.getConstrs():
            if c.IISConstr:
               print('%s' % c.constrName)
    end_time = time.time()
    elapsed_time = end_time - start_time
    print("Elapsed time: ", elapsed_time)

C = 6 #Number of customers/requests
Num_F = 2 #Number of SDLs
F = [i for i in range(C+1,C+1+Num_F)] #Set of SDL-Nodes
I = [i for i in range(1,C+1)] #Set of AHD-customer nodes
N = I + F #Set of AHD-customer nodes and SDL nodes
N_D = N.copy()# Set N with depot node
N_D.insert(0,0)
V = [] #Set of SDL-nodes for each request i
for v in I:
    F_cus_no = int(random.uniform(1,3))
    F_cus = random.sample(F,F_cus_no)
    F_cus.sort()
    V.append(F_cus)
delta = 1 #Compensation for delivery to SDL
gamma = 3 #Fixed Cost for each Vehicle used
B = [] #Capacity of SDL-Boxes
for b in F:
    B.append(int(random.uniform(1,3)))
#Distance matrix cij    
c = np.zeros((len(N_D),len(N_D)))
for i in N_D:
    for j in N_D:
        if j!=i:
            c[i,j] = int(random.uniform(5,8)) 
        else:
            c[i,j] = float("inf")
#travel time matrix tij
t = c #equals distance from i to j
e = [0]
e_c = [3 for i in I] #Service time of each customer i
e_f = [5 for f in F] #Service time of each SDL
e = e + e_c + e_f #Service time of customers and SDLs
Tmax = 35 #Maximum travel time for vehicles used 
a = [0 for i in I] #lower bound of time windows of customers i
a.insert(0,0) #earliest time where vehicles are allowed to leave the depot
l = [int(random.uniform(20,30)) for i in I] #upper bound of time windows of customers i 
l.insert(0,Tmax) #latest time to get back to the depot

Rundreise = VRPTW_SDL(N,N_D,I,F,V,e,Tmax,a,l,delta,gamma,t,c,C) 

The result:

    Variable            X 
-------------------------
      x[0,0]            1 
      x[0,5]            1 
      x[1,1]            1 
      x[2,2]            1 
      x[3,3]            1 
      x[4,4]            1 
      x[5,0]            1 
      x[5,5]            1 
      x[6,6]            1 
      Y[7,7]            1 
      Y[7,8]            1 
      Y[8,8]            1 
        s[5]            5 

14.0

I would be very happy if you could help me to find the reason for these wrong results.

Kind regards

Alex