OR-TOOLS : delivery node with multiple possible pickup nodes

I am using ortools to model a VRP with pickup and delivery constraints, where pickups can be done at different nodes. For example, if node A has a demand, it must be picked at node B or C.

Here is how I do this:

# data["pickups_deliveries"] is a dict with keys delivery_nodes and values a list of possible pickup nodes
# example : data["pickups_deliveries"][a] = [b,c]

for delivery_node in self.data["pickups_deliveries"]:
# choose one node among all pickup options
all_pickups = [
self.manager.NodeToIndex(p)
for p in self.data["pickups_deliveries"][delivery_node]
]

# same vehicle for pickup and delivery
delivery_index = self.manager.NodeToIndex(delivery_node)
sum(
self.routing.ActiveVar(p) * self.routing.VehicleVar(p)
for p in all_pickups
)
== self.routing.VehicleVar(delivery_index)
)

# precedence constraint
time_dimension = self.routing.GetDimensionOrDie("Time")
sum(
self.routing.ActiveVar(p)
* time_dimension.CumulVar(p)
for p in all_pickups
)
<= time_dimension.CumulVar(delivery_index)
)


This works if everything fits into 1 vehicle. But if capacity constraints require more than 1 vehicle, the solver does not find a solution after a few minutes (solver status 3).

I suspect something is wrong with the constraints imposing that the same vehicle is used for the pickup and delivery, but I am not sure. I have also tried using the following code (from here), but it does not help:

pickup_vehicles = [self.routing.VehicleVar(i) for i in all_pickups]
deliver_vehicle = [self.routing.VehicleVar(delivery_index)]