Apparent VRP Paradox

[juanluisrto]

I have run an experiment with pyvroom and have obtained some strange results.
I'd like to share them here in case someone knows why this happens or in case it opens the door to identifying some bug.

My experiment consists in solving a VRP in two different ways in order to compare their total cost / duration.

• I define a set of jobs at random locations with random weights and a set of vehicles also initialized at random (coordinates are sampled uniformly inside a circle)
• For the first way, I solve a regular VRP which we call simple_vrp = VRP(jobs, vehicles)
• For the second way:
  • I split my fleet of vehicles in two by half. So now I have vehicles_first, vehicles_second
  • I solve a VRP with all the jobs but only the first half of vehicles split_first_vrp = VRP(jobs, vehicles_first)
  • I retrieve the unassigned jobs from that VRP jobs_left = split_first_vrp.unassigned
  • I solve a second VRP with the other half of the fleet but only for the remaining jobs: split_second_vrp = VRP(jobs_left, vehicles_second)

So I am comparing the total cost of the simple_vrp vs the sums of the costs of split_first_vrp and split_second_vrp.
The paradox is that sometimes cost(simple_vrp) > cost(split_first_vrp) + cost(split_second_vrp)
This is completely counterintuitive to me, since the solver should be able to optimize the fleet usage when knowing all the vehicles available, compared to solving with a subset of the vehicles and then completing the jobs left with the rest of vehicles.

To reproduce this behaviour, here is some code:

import vroom
import numpy as np
import math
from scipy.spatial.distance import cdist


def generate_random_points( center : tuple, n: int, r: float, seed: int = 42):
    "Generates n uniformly distributed random points around a center point within a radius r"
    rng = np.random.default_rng(seed=seed)
    points = []
    cx, cy = center

    for _ in range(n):
        angle = rng.uniform(0, 2 * math.pi)
        distance = rng.uniform(0, r)
        dx = distance * math.cos(angle)
        dy = distance * math.sin(angle)
        new_point = (cx + dx, cy + dy)
        points.append(new_point)

    return points

N_JOBS = 400
N_VEHICLES = 8

SPLIT = int(N_VEHICLES/2) # index at which the fleet is split (half)

center = (0,0)

jobs_loc = generate_random_points(center, N_JOBS, r = 1)
vehicles_loc = generate_random_points(center, N_VEHICLES, r = 0.4)
points = jobs_loc + vehicles_loc

# Calculate the Euclidean distance matrix
fake_duration_matrix = cdist(points, points, metric='euclidean') * 100



def define_and_solve_vrp(jobs, vehicles):
    problem = vroom.Input()
    
    problem.set_durations_matrix(
        profile="car",
        matrix_input=fake_duration_matrix,
    )
    problem.add_job(jobs)
    problem.add_vehicle(vehicles)
    return problem.solve(exploration_level=5, nb_threads=4)
    

jobs = [vroom.Job(i, location = points.index(loc)) for i, loc in enumerate(jobs_loc)]
vehicles = [vroom.Vehicle(i,
                      start=points.index(loc),
                      end=points.index(loc)            
                     )
            for i, loc in enumerate(vehicles_loc)]
           
# We solve the VRP for all jobs and vehicles
simple_vrp_solution = define_and_solve_vrp(jobs, vehicles)

# We split the fleet
vehicles_first, vehicles_second = vehicles[:SPLIT], vehicles[SPLIT:]

# We solve the first part of the VRP with 
split_first_vrp_solution = define_and_solve_vrp(jobs, vehicles_first)

jobs_left = []

# We recreate the remaining jobs
for job in split_first_vrp_solution.unassigned:
    i = job.index()
    loc = jobs_loc[i]
    jobs_left.append(
        vroom.Job(i, location = points.index(loc))
    )

# We solve the second part of the VRP with the jobs left from the first part 
split_second_vrp_solution = define_and_solve_vrp(jobs_left, vehicles_second)


normal_vrp_cost = simple_vrp_solution.summary.cost 
split_vrp_cost = (split_first_vrp_solution.summary.cost + split_second_vrp_solution.summary.cost)

n_done_jobs_simple = N_JOBS - len(simple_vrp_solution.unassigned)
n_done_jobs_split  = N_JOBS - len(split_second_vrp_solution.unassigned)


print(f"normal_vrp_cost = {normal_vrp_cost} - Nº completed jobs = {n_done_jobs_simple}")
print(f" split_vrp_cost = {split_vrp_cost} - Nº completed jobs = {n_done_jobs_split}")
print(f"Paradox? = {normal_vrp_cost > split_vrp_cost}")

With this configuration the output of the program is:

>>>normal_vrp_cost = 2350 - Nº completed jobs = 400
>>> split_vrp_cost = 2340 - Nº completed jobs = 400
>>>Paradox? = True

There are many other configurations in which this happens (N_JOBS = 100, N_VEHICLES = 3), although it's not the usual case.
Do you know why this could be happening? Surely not related to pyvroom per se and its more on the vroom side

[jonathf]

I am not sure if I agree. Vroom is designed to optimize on the problem here and now, not on some potential larger question it doesn't know about yet. I'm not surpriced at all that sometimes the cost is a little higher like that.

That said, you are right that this is a vroom issue. @jcoupey do you want to add your two cents?

[jcoupey]

First, in principle there is no guarantee that you'll get normal_vrp_cost <= split_vrp_cost. This could indeed sound weird since the "normal" instance does have a more global view of the problem, but you'd only get this kind of assurance if we were to provide guaranteed optimal solutions.

That said, I suspect there is something else going on here related to the specific setup you're trying.

I define a set of jobs at random locations with random weights

Do you mean capacity restriction by "weights"? In your code there is no constraint whatsoever so basically there is nothing preventing any vehicle to handle all jobs in a very long route. Thus I'd be surprised if jobs_left is not empty for your second solving round.

Surely not related to pyvroom per se and its more on the vroom side

Yes. Is there a convenient way to serialize the problem instances to the "usual" json format from pyvroom? That would make things easier to actually look at the instances and reproduce.

[jonathf]

Pyvroom supports reading json format on input, and writing json for solutions. It does not support writing json for input. Is that supported in upstream vroom?

[juanluisrto]

Thanks for the insights to both of you, I understand that there are no guarantees of getting more optimal solutions with more vehicles. It just seemed to me weird enough to flag it.
Regarding the capacity restrictions, forget what I said. I removed those constraints and realized this behaviour still occurred (also removed max_travel_time), so now I'm just running unconstrained VRPs.

@jcoupey Its true that jobs_left was many times empty, unless the jobs were sufficiently close to the initial location of another vehicle.

I have now simplified my experiment and run some tests.
The new setup is the following, given a set of jobs and a set of vehicles:

1. Generate all possible fleets from the set of vehicles (any subset of 1 or more vehicles)
2. For each fleet, solve the VRP for all jobs, record its cost
3. identify the minimum cost to solve the VRP for 1,2,3 ... vehicles

I run this experiment 5 times with different amounts of jobs (N_JOBS = [10, 50, 100, 150, 200]) and fleets of up to 8 vehicles and here are the results.

Highlighted in yellow is you can see the minimum cost needed to solve each VRP (columns). It's notable that for the VRP with 150 Jobs, the solution with 1 vehicle outperforms any other solution with up to 8 vehicles.

[jcoupey]

Its true that jobs_left was many times empty, unless the jobs were sufficiently close to the initial location of another vehicle

I don't get that: if you have no constraint, any number of vehicles should be able to handle all jobs. If you have no constraint and jobs_left is not always empty then there is something else screwed in your setup.

It's notable that for the VRP with 150 Jobs, the solution with 1 vehicle outperforms any other solution with up to 8 vehicles.

Hard to draw any conclusion on optimization performance without knowing how your vehicles are defined:

• if all vehicles have the same start/end locations, then it's always cheaper to use a single vehicle if the constraints allows it;
• if vehicles are geographically distributed, then using more vehicles can lower the cost.

In the latter case, it's not a matter of having a "better" solution, it's simply about using vehicles that are closer to the jobs.

[jcoupey]

@jonathf no we don't have an ability to serialize the problem to json upstream because usually it's already defined as json in input (that is except for the folks using the C++ libvroom API, kind of a similar situation as the one with pyvroom).

[juanluisrto]

Thanks again @jcoupey for your answers.

I don't get that: if you have no constraint, any number of vehicles should be able to handle all jobs. If you have no constraint and jobs_left is not always empty then there is something else screwed in your setup.

I can confirm that jobs_left is always empty once I removed the max_travel_time and capacity constraints from the vehicles.

If vehicles are geographically distributed, then using more vehicles can lower the cost.

That is the case, I'm creating vehicles at random locations, with the default operating costs (zero fixed cost, only cost per hour). I still believe it is strange is that with a fraction of a fleet I get lower costs than with the whole fleet.

Here is the code for my simplified experiment.

import vroom
import numpy as np
import pandas as pd
import math
from scipy.spatial.distance import cdist

from tqdm import tqdm
import itertools


def generate_random_points( center : tuple, n: int, r: float, seed: int = 42):
    "Generates n uniformly distributed random points around a center point within a radius r"
    rng = np.random.default_rng(seed=seed)
    points = []
    cx, cy = center

    for _ in range(n):
        angle = rng.uniform(0, 2 * math.pi)
        distance = rng.uniform(0, r)
        dx = distance * math.cos(angle)
        dy = distance * math.sin(angle)
        new_point = (cx + dx, cy + dy)
        points.append(new_point)

    return points


def create_random_jobs_and_vehicles(N_JOBS, N_VEHICLES, seed = 42):
    center = (0,0)
    jobs_loc = generate_random_points(center, N_JOBS, r = 1)
    vehicles_loc = generate_random_points(center, N_VEHICLES, r = 0.4)
    points = jobs_loc + vehicles_loc

    duration_matrix = cdist(points, points, metric='euclidean') * 100

    jobs = [vroom.Job(i, location = points.index(loc)) for i, loc in enumerate(jobs_loc)]
    vehicles = [vroom.Vehicle(i,
                          start=points.index(loc),
                          end=points.index(loc),
                         )
                for i, loc in enumerate(vehicles_loc)]

    return jobs, vehicles, duration_matrix


def solve_vrp(jobs, vehicles, matrix):
    problem = vroom.Input()
    
    problem.set_durations_matrix(
        profile="car",
        matrix_input=matrix,
    )
    problem.add_job(jobs)
    problem.add_vehicle(vehicles)
    return problem.solve(exploration_level=5, nb_threads=4)


def vrp_with_all_fleet_combinations(jobs, vehicles, matrix):
    """Solves the same VRP for each possible subselection of the vehicles list.
    If vehicles = [A, B], it solves the VRP for [A], [B], [A,B]"""
    metrics = []
    for i in tqdm(range(1, len(vehicles) + 1), desc='Outer Loop'):
        for fleet in itertools.combinations(vehicles, i):
            id = "_".join([str(v.id) for v in fleet])
            solution = solve_vrp(jobs, fleet, matrix)
            metrics.append(
                {
                    "id": id,
                    "n_vehicles": len(fleet),
                    "n_jobs" : len(jobs),
                    "cost": solution.summary.cost,
                    
                }
            )
    return metrics


all_metrics = []

N_JOBS_LIST = [10, 50, 100, 150]
max_n_vehicles = 8

for n_jobs in N_JOBS_LIST:
    jobs, vehicles, duration_matrix = create_random_jobs_and_vehicles(n_jobs, max_n_vehicles, seed = 42)
    print(n_jobs, max_n_vehicles)
    metrics = vrp_with_all_fleet_combinations(jobs, vehicles, duration_matrix)
    all_metrics.extend(metrics)



df = (pd.DataFrame(all_metrics)
      .set_index("id")
      .groupby(["n_jobs", "n_vehicles"])
      .min()
      .reset_index()
      .pivot(columns = "n_jobs", values = "cost", index = "n_vehicles")
     )

df.style.background_gradient(axis = 0).highlight_min()

[juanluisrto]

Hi @jcoupey @jonathf , I would like to retake this conversation and ask for your input if its ok.

I have crafted a very simple setup in which the unexpected behaviour shows up.
I'd love to know a bit more of why this happens, either if it is to understand why my assumptions are wrong or if it helps to identify a possible improvement of vroom.

The setup is the following, I have 5 jobs and 2 vehicles. Vehicles should start and end in the same place. There is no service time, no time windows and no weight constraints.

First I solve the VRP with vehicle 1, then with vehicle 2, and then with both vehicles. The first VRP has a total duration of 388, whereas the others have a duration of 396. I expected the VRP with both vehicles to only use vehicle 1 and thus have also a cost of 388.

import vroom
import numpy as np

duration_matrix = \
np.array( [[  0.,  82.,  54., 135., 130.,  46.,  68.],
           [ 82.,   0.,  95., 134.,  91.,  73.,  50.],
           [ 54.,  95.,   0.,  84.,  97.,  22.,  52.],
           [135., 134.,  84.,   0.,  63.,  90.,  88.],
           [130.,  91.,  97.,  63.,   0.,  88.,  63.],
           [ 46.,  73.,  22.,  90.,  88.,   0.,  33.],
           [ 68.,  50.,  52.,  88.,  63.,  33.,   0.]])


jobs = [vroom.Job(i, location = vroom.LocationIndex(i)) for i in range(0,5)]

vehicle_1 = vroom.Vehicle(1, start=5, end=5)
vehicle_2 = vroom.Vehicle(2, start=6, end=6)



def solve_vrp(jobs, vehicles, matrix):
    problem = vroom.Input()
    
    problem.set_durations_matrix(
        profile="car",
        matrix_input=matrix,
    )
    problem.add_job(jobs)
    problem.add_vehicle(vehicles)
    return problem.solve(exploration_level=5, nb_threads=4)
            

solution_vehicle_1 = solve_vrp(jobs, [vehicle_1], duration_matrix)

solution_vehicle_2 = solve_vrp(jobs, [vehicle_2], duration_matrix)

solution_both_vehicles = solve_vrp(jobs, [vehicle_1, vehicle_2], duration_matrix)


print(f"duration with v1      = {solution_vehicle_1.summary.duration}")
print(f"duration with v2      = {solution_vehicle_2.summary.duration}")
print(f"duration with v1 & v2 = {solution_both_vehicles.summary.duration}")

Output
>>>duration with v1      = 388
>>>duration with v2      = 396
>>>duration with v1 & v2 = 396

In the third problem, when vroom has both vehicle_1 and vehicle_2 available to solve the VRP, it decides to assign all jobs to vehicle_2 instead of using vehicle_1 which would yield a more optimal route.

Is there something I am missing?

[jcoupey]

@juanluisrto really appreciate the efforts to narrow this down to a minimal example! The one you have here looks like it you be pretty straightforward to translate to the upstream json format. Do you think you could post that in a dedicated upstream ticket? This would both ease looking the problem up (at least for me), and allow to close this ticket that is not related to pyvroom per say any more.