VRPTW Infeasible

I was trying to code this model in python with gurobi, but I have some problems that caused the infeasible. Here is my code: 

from gurobipy import Model, GRB, quicksum
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import random
import math
import time

#generate data------------------------------------------------------------------------------
np.random.seed(5)

# Parameters
N = 5  # Number of regular stops
M = 2  # Number of landfills
K = 2  # Number of vehicles
C = 100  # Vehicle capacity

# Generating service times (s_i) for regular stops and landfills (randomly between 0 to 10), depot, and lunch break with 0
service_times = np.concatenate(([0], np.random.randint(0, 11, N + M), [0]))
# Earliest possible departure time from the depot
E = 0
# Latest possible arrival time at the depot
L = 600
# Generating time windows ([a_i, b_i]) for nodes, with the lunch break specifically between 12:00 and 14:00
time_windows_start = np.concatenate(([0], np.random.randint(1, 240, N + M), [4*60]))  # Start times
time_windows_end = np.concatenate(([600], np.random.randint(360, 600, N + M), [6*60]))  # End times

# Generating demands (d_i) for regular stops, with zero demand for landfills, lunch break, and depot
demands = np.concatenate(([0], np.random.randint(10, 30, N), np.zeros(M + 1)))

# Generating travel times (t_ij) between all pairs of nodes, symmetric and random between 1 to 20, zero for self-travel
travel_times = np.random.randint(1, 20, size=(N + M + 2, N + M + 2))
np.fill_diagonal(travel_times, 0)

# Compile generated data into a DataFrame for display
df_data = pd.DataFrame({
    "Service Times": service_times,
    "Time Window Start": time_windows_start,
    "Time Window End": time_windows_end,
    "Demands": demands
})

# Display the DataFrame to verify the generated data
print(df_data)

# Optionally, display the travel times matrix
print("Travel Times Matrix:")
print(travel_times)

#Preparing model------------------------------------------------------------------------------
model = Model("VRPTW")


# Decision Variables
# X[i, j, k]: 1 if vehicle k travels from node i to node j, 0 otherwise
X = model.addVars(N + M + 2, N + M + 2, K, vtype=GRB.BINARY, name="X")

# W[i, k]: Start service time at node i for vehicle k
W = model.addVars(N + M + 2, K, vtype=GRB.CONTINUOUS, name="W")

# D[i, k]: Cumulative demand at node i for vehicle k
D = model.addVars(N + M + 2, K, vtype=GRB.CONTINUOUS, name="D")

for i in range(N + M + 2):  
    for k in range(K):  
        D[i, k].lb = 0  
        D[i, k].ub = C  

# Objective Function: Minimize total travel time
model.setObjective(quicksum(X[i, j, k] * travel_times[i, j] 
                            for i in range(N + M + 2) 
                            for j in range(N + M + 2) 
                            for k in range(K)), GRB.MINIMIZE)

# Constraint (1): Each regular stop is served exactly once
for i in range(1, N + 1):
    model.addConstr(quicksum(X[i, j, k] for k in range(K) for j in range(N + M + 2) if j != i) == 1, f"visit_once_{i}")
    
# Constraint (2): Each vehicle starts from the depot
for k in range(K):
    model.addConstr(quicksum(X[0, j, k] for j in range(1, N + 1)) == 1, f"start_from_depot_{k}")
    
# Constraint (3): The cumulative demand at each stop for each vehicle must not exceed vehicle capacity
for k in range(K):
    for i in range(1, N + M + 1):  # Assuming that the depot and lunch break do not contribute to the load
        model.addConstr(D[i, k] <= C, f"capacity_{i}_{k}")
    
# Constraint (4): Cumulative demand is reset to zero at the depot and each landfill
for k in range(K):
    # Reset at the depot (node 0)
    model.addConstr(D[0, k] == 0, f"depot_reset_{k}")
    # Reset at each landfill node
    for m in range(N + 1, N + M + 1):
        model.addConstr(D[m, k] == 0, f"landfill_reset_{m}_{k}")
    
bigM = 1000  

# Constraint (5): Update cumulative demand correctly when a stop is visited
for k in range(K):
    for i in range(N + M + 2):
        for j in range(N + M + 2):
            if i != j:  
                model.addConstr(D[i, k] + demands[j] - D[j, k] <= (1 - X[i, j, k]) * bigM, f"cumulative_demand_{i}_{j}_{k}")

# Initialize Nk as integer decision variables for each vehicle
Nk = model.addVars(K, vtype=GRB.INTEGER, name="Nk")

# Landfill nodes indices
landfill_nodes = list(range(N + 1, N + M + 1))

# Add constraints to calculate Nk based on trips to landfill sites
for k in range(K):
    model.addConstr(
        Nk[k] == quicksum(X[i, j, k] for i in range(N + M + 2) for j in landfill_nodes),
        name=f"Nk_calculation_{k}"
    )

# Constraint (6): Ensure the total demand does not exceed the capacity times the number of disposal trips
for k in range(K):
    model.addConstr(
        quicksum(demands[j] * quicksum(X[i, j, k] for i in range(N + M + 2)) for j in range(1, N + 1)) 
        <= C * Nk[k], 
        f"capacity_trips_{k}"
    )

# Constraint (7): Number of trips from the regular stops and the depot to disposal facilities equals Nk
for k in range(K):
    model.addConstr(
        quicksum(quicksum(X[i, m, k] for m in range(N + 1, N + M + 1)) for i in range(N + M + 1)) 
        == Nk[k], 
        f"number_of_disposal_trips_{k}"
    )
    
# Constraint (8): The number of trips from disposal facilities to the regular stops and the depot is Nk - 1
for k in range(K):
    model.addConstr(
        quicksum(X[m, i, k] for m in range(N + 1, N + M + 1) for i in range(N + M + 2) if i != m) == Nk[k] - 1,
        f"trips_from_disposal_to_stops_{k}"
    )

# Constraint (9): The last trip from a disposal facility for each vehicle is to the depot
for k in range(K):
    model.addConstr(
        quicksum(X[m, 0, k] for m in range(N + 1, N + M + 1)) == 1,
        f"last_trip_to_depot_{k}"
    )    
    
# Constraint (10): Each vehicle arrives at the lunch break from one of the nodes
for k in range(K):
    model.addConstr(
        quicksum(X[i, N + M + 1, k] for i in range(1, N + M + 1)) == 1,
        f"arrive_at_lunch_{k}"
    )

# Constraint (11): Each vehicle departs from the lunch break to one of the nodes
for k in range(K):
    model.addConstr(
        quicksum(X[N + M + 1, j, k] for j in range(N + M + 1)) == 1,
        f"depart_from_lunch_{k}"
    )
    
# Constraint (12): A vehicle must leave a stop if it arrives at that stop
for k in range(K):
    for j in range(N + M + 2):  # Includes all nodes: stops, landfills, depot, and lunch break
        model.addConstr(
            quicksum(X[i, j, k] for i in range(N + M + 2)) == 
            quicksum(X[j, i, k] for i in range(N + M + 2)),
            f"flow_conservation_{j}_{k}"
        )
        
# Constraint (13): Time window and service time constraints without lunch break consideration
for k in range(K):
    for i in range(N + M + 2):
        for j in range(N + M + 2):
            if i != j:
                model.addConstr(W[i, k] + service_times[i] + travel_times[i, j] - W[j, k] <= (1 - X[i, j, k]) * bigM, f"time_window_{i}_{j}_{k}")

# Constraint (14): Time window and service time constraints with lunch break consideration
# This constraint will apply to the arcs leading into and out of the lunch break node
for k in range(K):
    for i in range(N + M + 2):
        for j in range(N + M + 2):
            if i != j:
                model.addConstr(W[i, k] + service_times[i] + service_times[N + M + 1] + travel_times[i, j] - W[j, k] <= (2 - X[i, N + M + 1, k] - X[N + M + 1, j, k]) * bigM, f"lunch_break_{i}_{j}_{k}")

# Constraint (15): Service start time must be within the time window of the node
for k in range(K):
    for i in range(N + M + 2):
        model.addConstr(time_windows_start[i] * quicksum(X[j, i, k] for j in range(N + M + 1)) <= W[i, k], f"service_start_earliest_{i}_{k}")
        model.addConstr(W[i, k] <= time_windows_end[i] * quicksum(X[i, j, k] for j in range(N + M + 1)), f"service_start_latest_{i}_{k}")

# Constraint (16): Service start time must be within the overall time horizon
for k in range(K):
    for i in range(N + M + 2):
        model.addConstr(E <= W[i, k], f"overall_horizon_start_{i}_{k}")
        model.addConstr(W[i, k] <= L, f"overall_horizon_end_{i}_{k}")
        
# Update the model with the new bounds
model.update()

#Optimization------------------------------------------------------------------------------
import time

# Optimize the model
start_time = time.time()  # Start the timer
model.optimize()  # Run the optimization
end_time = time.time()  # End the timer

# Calculate and print the running time
running_time = end_time - start_time
print(f"The optimization ran for {running_time:.4f} seconds")

# Check if the model has been solved to optimality
if model.status == GRB.OPTIMAL:
    print("Optimal solution found.")
    
    # Print the solution (for example, the service start times for each node and vehicle)
    for k in range(K):
        for i in range(N + M + 2):
            if W[i, k].X > 0:  # Only print times for nodes that are visited
                print(f"Vehicle {k} starts service at node {i} at time {W[i, k].X}")

    # print the routes by checking the values of X[i, j, k]
    for k in range(K):
        for i in range(N + M + 2):
            for j in range(N + M + 2):
                if X[i, j, k].X > 0.5:  # Assuming a binary variable, we check if it's 'active'
                    print(f"Vehicle {k} travels from node {i} to node {j}")
                    
elif model.status == GRB.INFEASIBLE:
    print("The model is infeasible; no solution exists.")
elif model.status == GRB.UNBOUNDED:
    print("The model is unbounded; no optimal solution can be found.")
else:
    print(f"Optimization was stopped with status {model.status}")