VRP with Simultaneaous Delivery and Stochastic Pickup with Timewindows

Hello, I am trying to implement VRPSPDTW with stochastic pickups using multi-objective modeling on Gurobi where one objective minimizes the total distance travelled and other onw maximizes the stochastic pickup based on probability of high, medium and low scenarios. For some reason, I am not able to integrate the BigM time constraint here, if I remove it, the model successfully does pickups, but has subtours, however after adding the BigM time constraint, the subtours are eliminated but no pickups are generated in the add. The code is as follows please:

import numpy as np
import gurobipy as gp
from gurobipy import GRB

# Random seed for reproducibility
np.random.seed(42)

m = gp.Model("VRPSPDTW")

# Number of nodes
N = 5

# Number of cars
K = 2

# Capacity of each car
C = 25

# Arcs
A = [(i,j) for i in range(N) for j in range(N) if i != j]

# Create xij index
x = m.addVars(A, vtype=GRB.BINARY, name='x')

d = {}
# Nodes coordinates and distances between them
Nodes = np.array([
    (35, 35),
    (41, 49),
    (35, 17),
    (55, 45),
    (55, 20)
])

for (i, j) in A:
    d[(i, j)] = np.linalg.norm([Nodes[i][0] - Nodes[j][0], Nodes[i][1] - Nodes[j][1]])

# Print the distance matrix
print("Distance Matrix:")
for i in range(N):
    for j in range(N):
        if (i, j) in d:
            print(f"d[{i},{j}] = {d[i, j]:.2f}", end='\t')
        else:
            print("d[{},{}] = N/A".format(i, j), end='\t')
    print()  # Newline for each row

# Delivery demand
quantity = np.array([0, 10, 7, 13, 19])

# Base Pickup demand
base_backhaul = np.array([0, 10, 7, 13, 19])

# Time windows (for simplicity, assume they are the same for all nodes and vehicles)
# a = np.array([0, 161, 50, 116, 149])  # earliest start time
# b = np.array([230, 171, 60, 126, 159])  # latest finish time

a = np.array([0, 10, 20, 30, 40])  # earliest start time
b = np.array([200, 120, 130, 140, 150])  # latest finish time

# Service times
service_times = np.array([0, 10, 10, 10, 10])

# Scenarios
scenarios = {
    'high': {'probability': 0.4, 'multiplier': np.random.uniform(1.2, 1.3)},
    'medium': {'probability': 0.3, 'multiplier': np.random.uniform(0.95, 1.05)},
    'low': {'probability': 0.3, 'multiplier': np.random.uniform(0.6, 0.8)}
}

# Create Rij index
R = m.addVars(A, vtype=GRB.INTEGER, name='R')

# Create Pij index
P = m.addVars(A, vtype=GRB.INTEGER, name='P')

# Create min_pickup index
min_pickup = m.addVars(N, scenarios.keys(), vtype=GRB.INTEGER, name='min_pickup')

# Create time variables
T = m.addVars(N, vtype=GRB.CONTINUOUS, name='T')

# Initialize Big M variable
M = {}
for (i, j) in A:
    M[i, j] = b[i] + service_times[i] + d[i, j] - a[j]

# Parameters
alpha = 20
beta = 800

# Objective 1: Maximize total collected quantity
objective_pickup = gp.LinExpr()
for scenario_name, scenario in scenarios.items():
    p_omega = scenario['probability']
    y_omega = scenario['multiplier']
    for i in range(1, N):
        for j in range(1, N):
            if i != j:
                objective_pickup += p_omega * y_omega * P[i, j]

m.setObjectiveN(-objective_pickup, index=0, priority=1, weight=beta, name="MaximizePickup")

# Objective 2: Minimize total traveled distance
objective_distance = sum(d[i, j] * x[i, j] for (i, j) in A)
m.setObjectiveN(objective_distance, index=1, priority=0, weight=alpha, name="MinimizeDistance")

# Add Constraints
# Summation of Xij = 1, where j is from 0 to N, for all i
m.addConstrs((sum(x[i, j] for i in range(N) if (i, j) in A) == 1 for j in range(1, N)), name='E1')

# Summation of Xij = 1, where i is from 0 to N, for all j
m.addConstrs((sum(x[i, j] for i in range(N) if (i, j) in A) == sum(x[j, i] for i in range(N) if (j, i) in A) for j in range(1, N)), name='E2')

# Summation Rij - qj = Summation Rji
m.addConstrs((sum(R[i, j] for i in range(N) if (i, j) in A) - quantity[j] == sum(R[j, i] for i in range(N) if (j, i) in A) for j in range(1, N)), name='E3')

# Constraints to define min_pickup
for scenario_name, scenario in scenarios.items():
    y_omega = scenario['multiplier']
    for j in range(1, N):
        m.addConstr(min_pickup[j, scenario_name] <= quantity[j])
        m.addConstr(min_pickup[j, scenario_name] <= base_backhaul[j] * y_omega)

# Flow conservation for pickups under each scenario
for scenario_name, scenario in scenarios.items():
    y_omega = scenario['multiplier']
    for j in range(1, N):
        m.addConstr(sum(P[i, j] * y_omega for i in range(N) if (i, j) in A) - min_pickup[j, scenario_name] == sum(P[j, i] * y_omega for i in range(N) if (j, i) in A), name='E4')

# Capacity constraints under each scenario
for scenario_name, scenario in scenarios.items():
    y_omega = scenario['multiplier']
    m.addConstrs((R[i, j] + P[i, j] * y_omega <= C * x[i, j] for i in range(N) for j in range(N) if i != j), name=f'E7_{scenario_name}')

# Summation Ri0 = 0
m.addConstrs(((sum(R[i, 0] for i in range(1, N) if (i, 0) in A) == 0) for j in [0]), name='E5')

# Summation P0j = 0
for scenario_name, scenario in scenarios.items():
    y_omega = scenario['multiplier']
    m.addConstrs(((sum(P[0, j] * y_omega for j in range(1, N) if (0, j) in A) == 0) for i in [0]), name='E6')

# Ensure that vehicles leave the depot and return
m.addConstrs((sum(x[0,i] for i in range(1,N)) == K for j in [0]), name='E10')
m.addConstr(sum(x[i, 0] for i in range(1, N) if (i, 0) in A) <= K, name='E11')

# Rij >= 0
m.addConstrs((R[i, j] >= 0 for i in range(N) for j in range(N) if i != j), name='E8')

# Pij >= 0
for scenario_name, scenario in scenarios.items():
    y_omega = scenario['multiplier']
    m.addConstrs((P[i, j] * y_omega >= 0 for i in range(N) for j in range(N) if i != j), name='E9')

# Time window constraints

m.addConstrs((T[i] + d[i, j] + service_times[i] <= T[j] + M[i, j] * (1 - x[i, j]) for (i, j) in A if i != j), name='E12') # by BigM

m.addConstrs((T[i] >= a[i] for i in range(N)), name='E13_lower')
m.addConstrs((T[i] <= b[i] for i in range(N)), name='E13_upper')

# Optimize model
m.optimize()

# Check for infeasibility and print IIS if found
if m.status == GRB.INFEASIBLE:
    m.computeIIS()
    print('\nThe following constraints and variables are in the IIS:')
    for c in m.getConstrs():
        if c.IISConstr:
            print('%s' % c.constrName)
else:
    # Print decision variables and objective values
    for v in m.getVars():
        if v.X != 0:
            print(f'{v.varName}: {v.X}')

    # Print objective values
    for i in range(m.NumObj):
        print(f'Objective {i}: {m.getObjective(i).getValue()}')

    # Calculate and print total Rij and Pij fulfilled
    total_R = sum(R[i, j].X for i, j in A if R[i, j].X > 0)
    total_P = sum(P[i, j].X for i, j in A if P[i, j].X > 0)
    print(f'Total R (Deliveries) fulfilled: {total_R}')
    print(f'Total P (Pickups) fulfilled: {total_P}')

# Retrieve the optimal route
optimal_route = [(i,j) for (i,j) in A if x[i,j].x > 0.5]

# Sort the optimal route
start_location = 0  # assuming the route starts at location 0
sorted_route = []
while optimal_route:
    for i, route in enumerate(optimal_route):
        if route[0] == start_location:
            start_location = route[1]
            sorted_route.append(route)
            optimal_route.pop(i)
            break

# Create an array of locations in the sorted route
route_array = [route[0] for route in sorted_route]
route_array.append(sorted_route[-1][1])  # add the final location
route_array.pop(0)
print(route_array)

# Split the route array into separate routes
routes = []
current_route = []
for i, location in enumerate(route_array):
    current_route.append(location)
    if location == start_location and i != 0:
        routes.append(current_route)
        current_route = []

# Add the last route
if current_route:
    routes.append(current_route)

# Prepend the start location to each route
for i, route in enumerate(routes):
    if i > 0:
        route.insert(0, start_location)

# Print the routes
for i, route in enumerate(routes):
    print(f"Route {i+1}:")
    print(route)

The screenshot of mathematical model is also attached. Can someone please help me here urgently? I have been trying to fix this since quite many days with no luck... :(

One correction: beta is far greater than alpha since we want to prioritize pickup collection.