Prize Collecting VRP returns empty solution

[rastislavkopal]

What is the correct usage?

I am implementing a Prize Collecting VRP, where I have list of clients with their geoLocation. Each client has a prize that we will collect, and the service time (duration spent at a client). I have one vehicle that needs to maximize the prize collected within a day. Each client will also have a time window (availability time). We need to maximize the prize collected.

# Clients: name, coordinates, prize, duration (service time)
CLIENTS_DATA = [
    ["Client 1", (48.142, 17.099), 10, 120],
    ["Client 2", (48.1445, 17.1074), 8, 60],
    ["Client 3", (48.1429, 17.1133), 7, 30],
    ["Client 4", (48.1736, 16.9825), 9, 90],
    ["Client 5", (48.1447, 17.1075), 7, 60],
    ["Client 6", (48.1405, 17.1072), 8, 120],
    ["Client 7", (48.1348, 17.1015), 5, 90],
]

as the library does not support float coordinates, I tried to just multiply location by a constant (1000) to get INTs. I calculated the distances with geopy library. Now I am using UTM representation instead and just calculate the Manhattan distance between them. I actually do not know what is the standard way to encode the geoLocations as INT the best way.

Now when I create the model with the scores, coordinates, there is no found solution

def convert_latlon_to_utm(lat, lon):
    easting, northing, zone_number, zone_letter = utm.from_latlon(lat, lon)
    return int(easting), int(northing) # return in meters

START_INDEX = 0
COORDS = [(convert_latlon_to_utm(client[1][0], client[1][1])) for client in CLIENTS_DATA]
PRIZES = [client[2] for client in CLIENTS_DATA]
PRIZES[START_INDEX] = 0  # starting point has no prize
TIME_WINDOWS = [(0, 24 * 60) for _ in range(len(CLIENTS_DATA))] # for now, time window is the whole day

The model:

m = Model()

# We have only one depot with opening at 8:00 and closing at 16:00 
# What does early and late mean in the context od depot? That the vehicles can start 
# at the tw_early and must return before tw_late?
depots = [
    m.add_depot(
        x=COORDS[START_INDEX][0],
        y=COORDS[START_INDEX][1],
        tw_early= 8 * 60,
        tw_late =16 * 60,
    )
]

m.add_vehicle_type(1, depot=depots[0], max_duration=10000)

# Now we set the time windows, prizes and coordinates
clients = [
    m.add_client(
        x=COORDS[0][0],
        y=COORDS[0][1],
        tw_early=TIME_WINDOWS[0][0],
        tw_late=TIME_WINDOWS[0][1],
        prize=PRIZES[0],
        required=False,
    )
]

locations = depots + clients
for frm_idx, frm in enumerate(locations):
    for to_idx, to in enumerate(locations):
        distance = abs(frm.x - to.x) + abs(frm.y - to.y)  # Manhattan
        m.add_edge(frm, to, distance=distance)

res = m.solve(stop=MaxRuntime(5)) 

The solver returns either an empty solution or sometimes the warning:

An empty solution is being added to the population. This typically indicates that there is a significant difference between the values of the prizes and the other objective terms, which hints at a data problem. Note that not every part of PyVRP can work gracefully with empty solutions.

What can be the reason ? Thank you for any insights.

[N-Wouda]

Hi @rastislavkopal! 👋

Your prize vector has values between 5 and 10. Your distances, in meters, are orders of magnitude larger. PyVRP optimizes an objective of the form "minize distance - prizes" (the docs contain the exact formulation). To get a prize of 5 to 10, a vehicle needs to spend many thousands in distance to get there. That's not a good decision to make, so the algorithm decides that not visiting anything is the best it can do.

The solution here is to also scale the prizes, to make their values comparable to the distance values. That should make the prizes interesting enough for PyVRP to consider visiting the clients.

[rastislavkopal]

Oh, I see. Thanks for pointing that out @N-Wouda .

Now when I scale the prizes (I tried multiple scaling_factors with no success)

def scale_scores(scores, distance_matrix, scale_factor=1):
    distances = [distance for row in distance_matrix for distance in row if distance != 0] # Flatten and remove 0 

    median_distance = sorted(distances)[len(distances) // 2]
    median_price = sorted(scores)[len(scores) // 2]

    scale_factor = (scale_factor * median_distance) / median_price
    return [round(prize * scale_factor) for prize in scores]

The solution visits everything if I do not set any duration to edges and max_duration to vehicle. That makes sense.

When I set both duration and distance to each edge:

for frm_idx, frm in enumerate(locations):
    for to_idx, to in enumerate(locations):
        if frm_idx == to_idx:
            continue
        distance = DISTANCE_MATRIX[frm_idx][to_idx]
        walking_duration = round(distance / 83.33) 
        m.add_edge(frm, to, distance=distance, duration=walking_duration)

And max distance to the only vehicle I have:

m.add_vehicle_type(1, depot=depots[0], max_duration=12 * 60) 

There is always just INFEASIBLE solutions. Am I still missing something?

[N-Wouda]

@rastislavkopal could you provide me with a full script that I can run locally? I find it somewhat difficult to piece this together into a full, running example. Just copy pasting your local script should be sufficient, if it is not too large.

[N-Wouda]

@rastislavkopal you're running into an issue that's been reported previously in #461 and #508. Basically, our default time windows are [0, 0] for depots and clients. That's been fixed in main to [0, +inf), and will be part of 0.8.0 when we release that version in a few weeks.

For now, if you want to use durations, I suggest you provide time windows for depots and clients explicitly. In your sample something like tw_early=0 and tw_late=24 * 60 should work. With those settings I find a feasible solution within a second:

```python from pyvrp import Model from pyvrp.stop import MaxRuntime import numpy as np from geopy.distance import geodesic LANDMARKS = [ ["Statue of Liberty", (40.6892, -74.0445), 10, 60], ["Central Park", (40.785091, -73.968285), 4, 60], ["Times Square", (40.7580, -73.9855), 3, 30], ["Empire State Building", (40.748817, -73.985428), 5, 90], ["Brooklyn Bridge", (40.7061, -73.9969), 6, 30], ["Metropolitan Museum of Art", (40.779437, -73.963244), 1, 120], ["One World Observatory", (40.712743, -74.013379), 10, 90], ["High Line", (40.7480, -74.0048), 4, 45], ["9/11 Memorial & Museum", (40.7115, -74.0134), 2, 60], ["Broadway Theater District", (40.7580, -73.9855), 1, 120], ["The Museum of Modern Art", (40.7614, -73.9776), 8, 90], ["The Frick Collection", (40.7712, -73.9673), 4, 60], ["New York Public Library", (40.7532, -73.9822), 7, 45], ["Bryant Park", (40.7536, -73.9832), 3, 60], ["American Museum of Natural History", (40.7813, -73.9739), 9, 150], ["Greenwich Village", (40.7336, -74.0027), 5, 120], ["Chelsea Market", (40.7420, -74.0048), 4, 90], ["Battery Park", (40.7033, -74.0170), 2, 60], ["Wall Street", (40.7074, -74.0113), 3, 30], ] def calculate_distance_matrix(landmarks, unit='m', round_to=0): num_landmarks = len(landmarks) distance_matrix = [[0 for _ in range(num_landmarks)] for _ in range(num_landmarks)] for i in range(num_landmarks): for j in range(num_landmarks): if i != j: if unit == 'm': distance = round(geodesic(landmarks[i][1], landmarks[j][1]).meters, round_to) distance_matrix[i][j] = distance else: distance_matrix[i][j] = 0 return distance_matrix DISTANCE_MATRIX = calculate_distance_matrix(LANDMARKS, 'm', 2) ORIGINAL_PRIZES = [landmark[2] for landmark in LANDMARKS] # TODO balance the scale_factor def scale_scores(scores, distance_matrix, scale_factor=1): """Scale the scores of the landmarks by the distance from the starting point scale_factor: multiplier to enhance the median prize's influence in scaling """ distances = [distance for row in distance_matrix for distance in row if distance != 0] # Flatten the distance matrix and remove 0 # average_distance = sum(distances) / len(distances) median_distance = sorted(distances)[len(distances) // 2] # average_price = sum(scores) / len(scores) median_price = sorted(scores)[len(scores) // 2] scale_factor = (scale_factor * median_distance) / median_price return [round(prize * scale_factor) for prize in scores] PRIZES = scale_scores(ORIGINAL_PRIZES, DISTANCE_MATRIX) START_INDEX = 0 GEO_PRECISION = 1000 WALKING_SPEED_M_PER_MIN = 83.33 # Approx. 5 km/h in meters per minute DRIVING_SPEED_M_PER_MIN = 400 # Approx. 24 km/h in meters per minute, considering city driving TIME_SCALE = 60 m = Model() depots = [ m.add_depot( x=round(LANDMARKS[START_INDEX][1][0] * GEO_PRECISION), y=round(LANDMARKS[START_INDEX][1][1] * GEO_PRECISION), tw_early=0, tw_late=24 * TIME_SCALE, ) ] m.add_vehicle_type(1, depot=depots[0], max_duration=12*TIME_SCALE, name='Charlie') clients = [ m.add_client( x=round(LANDMARKS[idx][1][0] * GEO_PRECISION), y=round(LANDMARKS[idx][1][1] * GEO_PRECISION), tw_early=0, tw_late=24 * TIME_SCALE, prize=PRIZES[idx], required=False, ) for idx in range(1, len(LANDMARKS)) ] locations = depots + clients for frm_idx, frm in enumerate(locations): for to_idx, to in enumerate(locations): if frm_idx == to_idx: continue distance = DISTANCE_MATRIX[frm_idx][to_idx] duration = round(distance / WALKING_SPEED_M_PER_MIN) m.add_edge(frm, to, distance=distance, duration=duration) res = m.solve(stop=MaxRuntime(1)) print(res) ```

[N-Wouda]

@rastislavkopal does that resolve the issue?

[rastislavkopal]

Yes it did actually, thanks so much for that @N-Wouda .

[N-Wouda]

@rastislavkopal great!

[rastislavkopal]

I was also wondering if it is possible to add something like a vehicle breaks which must happen during the day. Let's say one hour for a lunch break within a specified time window.

For example using or-tools: https://github.com/google/or-tools/blob/stable/examples/cpp/cvrptw_with_breaks.cc

Btw. In your example in docs, you use

res = m.solve(stop=MaxRuntime(1), display=False) # one second

But it throws an error

TypeError: Model.solve() got an unexpected keyword argument 'display

[N-Wouda]

@rastislavkopal you're looking at the dev docs for 0.8.0, not the docs for the latest release (0.7.0). Those are here. We still need to update pyvrp.org to default to the latest stable release rather than the development version.

[rastislavkopal]

I really appreciate your help @N-Wouda . I know it's probably outside of the scope of the library and this discussion. But do you think it's possible to model the break for the driver during the day? As I mention, I need to set a for example 60 minute lunch break anywhere in time window 11:00 to 15:00. I was thinking about adding dummy nodes with zero distances which will be required. Not sure what is the standard way to handle this scenario.

[N-Wouda]

I was thinking about adding dummy nodes with zero distances which will be required.

This is probably your best bet for now. We do not have any formal support for lunch breaks (#352) or time-based breaks (#415) yet.