Vehicle Routing Problem with Fair Profits and Time Windows (VRP-FPTW)#

vrp_fptw.ipynb Open In Colab Open In Deepnote Open In Kaggle Open In Gradient Open In SageMaker Studio Lab Powered by AMPL

Description: This notebook implements and solves the Vehicle Routing Problem with Fair Profits and Time Windows (VRP-FPTW), a realistic and recent extension of the classical VRP problem.

A MILP formulation is used and implemented with Ampl, taking advantage of its automatic reformulation mechanism to write complex constraints like the time windows and the objective function to maximize the worst-off vehicule.

Tags: amplpy, routing, vrp, time-windows, mp, cuopt, highs, gurobi, research, automatic-reformulation, benchmark

Model author: Aitor Lopez <aitor.lopez@urjc.es>

Notebook author: Aitor Lopez <aitor.lopez@urjc.es> , Marcos Dominguez Velad <marcos@ampl.com>

References:

  • Sánchez, A. L., Lujak, M., Semet, F., & Billhardt, H. (2023, October). Vehicle routing problem with fair profits and time windows (vrp-fptw). In 2023 IEEE International Conference on Systems, Man, and Cybernetics (SMC) (pp. 1530-1536). IEEE.

  • Instances from: https://github.com/aitorls/DVRP-PFPTW-instances

🚚 1. What is VRP-FPTW?#

The Vehicle Routing Problem with Fair Profits and Time Windows (VRP-FPTW) is a new variant of the classic Vehicle Routing Problem (VRP), specifically designed for settings like crowdsourced delivery fleets or cooperative robot fleets.

Its goal is to:

  • Visit all customers exactly once, within their time windows.

  • Respect vehicle-specific constraints (capacity and autonomy).

  • Maximize the profit of the worst-performing vehicle in the fleet — promoting fairness across all participating vehicles.

This problem is particularly suitable for multi-agent environments, where:

  • Vehicles are independently managed (e.g., by different stakeholders).

  • Vehicles have private information (e.g., individual costs, speed, capacities).

  • Fair revenue-sharing is essential to encourage collaboration.

⚙️ 2. How does it work?#

The VRP-FPTW operates as follows:

✅ Inputs:#

  • Customers each with:

    • A demand to fulfill.

    • A revenue if served.

    • A time window \([l_i, u_i]\) during which they must be served.

  • Fleet of vehicles, each with:

    • A cost per distance traveled.

    • A maximum capacity (e.g., weight or volume).

    • An autonomy (e.g., max distance it can travel before needing recharge/refuel).

🎯 Objective:#

Instead of maximizing total profit or minimizing cost, the goal is to maximize the minimum profit among all vehicles (i.e., ensure the worst-off vehicle makes as much as possible).

Profit = Revenue from served customers - Cost of traveling.

3. MILP (Mixed Integer Linear Programming) Formulation#

The problem is mathematically formulated as a MILP model.

Parameters:#

Symbol

Description

\(N\)

Set of customers

\(V = \{0\} \cup N\)

Vertices including depot \(0\) and customers

\(K\)

Set of vehicles

\(d_{ij}\)

Distance between nodes \(i\) and \(j\)

\(c_{ijk}\)

Travel cost for vehicle \(k\) from \(i\) to \(j\)

\(t_{ijk}\)

Travel time for vehicle \(k\) from \(i\) to \(j\)

\(q_i\)

Demand of customer \(i\)

\(r_i\)

Revenue from serving customer \(i\)

\([l_i, u_i]\)

Time window for customer \(i\)

\(D_k\)

Autonomy of vehicle \(k\)

\(Q_k\)

Capacity of vehicle \(k\)

\(M\)

Large constant for time window constraints

Decision Variables:#

Variable

Meaning

\(x_{ijk}\)

1 if vehicle \(k\) travels from \(i\) to \(j\); 0 otherwise

\(v_i \geq 0\)

Time at which customer \(i\) is visited

\(y \in \mathbb{R}\)

Profit of the worst-off vehicle (to be maximized)

📏 Constraints:#

  1. Every customer is visited once:

\[ \sum_{k} \sum_{i} x_{ijk} = 1 \quad \forall j \in N \]

  1. Flow conservation for each vehicle:

\[ \sum_{i} x_{ijk} - \sum_{h} x_{jhk} = 0 \quad \forall j \in V,\ \forall k \in K \]

  1. Each vehicle used at most once:

\[ \sum_{j} x_{0jk} \leq 1 \quad \forall k \in K \]

  1. Respect vehicle autonomy:

\[ \sum_{i} \sum_{j} d_{ij} x_{ijk} \leq D_k \quad \forall k \in K \]

  1. Respect vehicle capacity:

\[ \sum_{i} \sum_{j} q_j x_{ijk} \leq Q_k \quad \forall k \in K \]

  1. Time window consistency:

Time window can be easily formulated via indicator constraints: $\( x_{i,j,k} = 1 \implies v_i + t_{ijk} \leq v_j \quad \forall (i,j),\ \forall k \)$

Alternative formulation using big-M: $\( v_i + t_{ijk} \leq v_j + M(1 - x_{ijk}) \quad \forall (i,j),\ \forall k \)$

  1. Time window satisfaction:

\[ l_i \leq v_i \leq u_i \quad \forall i \in N \]

  1. Variable domains:

\[ x_{ijk} \in \{0,1\},\quad v_i \in \mathbb{R},\quad y \in \mathbb{R} \]

Objective Function:#

We can define the objective through a non-linear expression using \(\min\) that’s automatically handle by Ampl:

\[\max \min_{k \in K} \sum_{i} \sum_{j} (r_i - c_{ijk}) x_{ijk} K\]

Maximize the minimum profit among all vehicles.

🔁 How It Is Solved:#

  • Centralized approach: Solves the full MILP directly, but becomes inefficient as the number of customers/vehicles increases.

  • Distributed MAS (Multi-Agent System):

    • Each vehicle is modeled as an autonomous agent solving a subproblem.

    • A fleet coordinator agent uses column generation to combine local solutions into a global one.

    • Vehicles don’t share private information (e.g., cost structures), only new profitable routes and marginal prices (shadow prices).

    • The system iteratively improves routes and assignments until convergence.

4. Implementation#

In this notebook, we are just solving the Centralized approach through a MILP formulation. This is not the most efficient way to solve the problem, but an illustrative case of how to formulate a VRP problem with Ampl in a solver-agnostic approach.

Different solvers can be used to solve the problem with open-source solvers like CuOpt from Nvidia or HiGHS, or commercial ones like Gurobi.

The MILP method is useful for small-size problems, research, or tools comparisons, but more advanced algorithms are explored in the original paper, and they should be used for larger instances.

We are using the public static instances for VRP-FPTW found in the repository:

https://github.com/aitorls/DVRP-PFPTW-instances

# Install dependencies
%pip install -q amplpy pandas numpy

# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook

ampl = ampl_notebook(
    modules=["highs", "cuopt", "gurobi"],  # modules to install
    license_uuid="default",  # license to use
)  # instantiate AMPL object and register magics

Model formulation with Ampl#

%%writefile vrp-fptw.mod
# Sets
set K;  # set of vehicles
set V;  # set of vertices (including depot)
set N within V;  # set of customers (N subset of V)
set A within V cross V;  # arcs (i,j)

# Parameters
param d{V, V};           # distance from i to j
param t{V, V};           # travel time from i to j
param c{V, V, K};        # cost of arc (i,j) for vehicle k
param r{V};              # revenue for visiting node j
param q{V};              # demand at node j
param Q{K};              # capacity of vehicle k
param D{K};              # distance limit for vehicle k
param l{V};              # lower time window at node
param u{V};              # upper time window at node

# Decision Variables
var x{V, V, K} binary;      # 1 if arc (i,j) used by vehicle k
var v{V};                   # arrival time at node j
var y;                      # objective value
var TotalProfit = sum{k in K, i in V, j in V} (r[i] - c[i,j,k]) * x[i,j,k]; # non-used

# Objective Function
maximize MinProfit:
    y;

# Constraints
s.t. Not_Enter_Same_Node{i in V, k in K}:
    x[i,i,k] = 0;

s.t. One_Visit{j in N}:
    sum{k in K, i in V} x[i,j,k] = 1;

s.t. Flow_Conservation{j in V, k in K}:
    sum{i in V} x[i,j,k] - sum{h in V} x[j,h,k] = 0;

s.t. Depot_Start{k in K}:
    sum{j in N} x[0,j,k] <= 1;

s.t. Distance_Limit{k in K}:
    sum{i in V, j in V} d[i,j] * x[i,j,k] <= D[k];

s.t. Capacity_Limit{k in K}:
    sum{i in V, j in N} q[j] * x[i,j,k] <= Q[k];

s.t. Time_Constraint{k in K, (i,j) in A: j != 0 and j in N}:
    x[i,j,k] == 1 ==> v[i] + t[i,j] <= v[j];

s.t. Time_Windows{i in V}:
    l[i] <= v[i] <= u[i];

s.t. Objective_Def:
    y = min{k in K} sum{i in V, j in V} (r[j] - c[i,j,k]) * x[i,j,k];

Overwriting vrp-fptw.mod

Prepare data#

import os
import requests
import numpy as np
import pandas as pd


def read_instance(filename):
    if os.path.exists(f"./instances/static/{filename}.npy"):
        filename = f"./instances/static/{filename}.npy"
    else:
        # Example: Downloading a notebook from GitHub (raw URL)
        instance_url = f"https://raw.githubusercontent.com/aitorls/DVRP-PFPTW-instances/main/instances/static/{filename}.npy"
        filename = f"{filename}.npy"

        response = requests.get(instance_url)
        if response.status_code == 200:
            with open(filename, "wb") as f:
                f.write(response.content)
            print(f"Downloaded: {filename}")
        else:
            print("Failed to download notebook.")

    data = np.load(filename, allow_pickle="TRUE").item()
    # Derive sets
    n_customers = data["n_customers"]
    V = list(range(n_customers + 1))
    N = list(range(1, n_customers + 1))
    K = list(range(data["vehicles"]))

    # Derive distances
    coords = data["node_coord"]
    edge_weights = np.linalg.norm(coords[:, None, :] - coords[None, :, :], axis=-1)

    # Prepare model parameter dictionaries
    d = data["edge_weight"]
    t = d.copy()
    c = {(i, j, k): edge_weights[i][j] for i in V for j in V for k in K}
    q = {j: int(data["demand"][j]) for j in V}
    r = {j: int(data["revenue"][j]) for j in V}
    l = {j: data["time_window"][j][0] for j in V}
    u = {j: data["time_window"][j][1] for j in V}
    Q = {k: data["capacity"] for k in K}
    D = {k: data["autonomy"] for k in K}
    A = [(i, j) for i in V for j in V if i != j]

    return n_customers, V, N, K, d, t, c, q, r, l, u, Q, D, A

Read the model and load data#

def read_model():
    model = AMPL()
    model.read("vrp-fptw.mod")
    return model


def data_to_ampl(model, instance):
    n_customers, V, N, K, d, t, c, q, r, l, u, Q, D, A = read_instance(
        filename=instance
    )

    # Load sets
    model.set["K"] = K
    model.set["V"] = V
    model.set["N"] = N
    model.set["A"] = A

    # Load parameters
    model.param["d"] = d
    model.param["t"] = t
    model.param["c"] = c
    model.param["q"] = q
    model.param["r"] = r
    model.param["l"] = l
    model.param["u"] = u
    model.param["Q"] = Q
    model.param["D"] = D

    # Set the value for big M
    # model.param['M'] = 1e6
    return model


def solve(instance, solver):
    model = read_model()
    data_to_ampl(model=model, instance=instance)
    # See solver options at https://dev.ampl.com/solvers/
    model.solve(solver=solver, mp_options="outlev=1 timelimit=150")
    return model


instance = "SFPTW_25_5_1"
# instance='SFPTW_100_20_0'
solver = "highs"

m = solve(instance=instance, solver=solver)
assert m.solve_result == "limit", m.solve_result

Downloaded: SFPTW_25_5_1.npy
HiGHS 1.10.0:   tech:outlev = 1
  lim:time = 150

AMPL MP initial flat model has 3277 variables (1585 integer, 1665 binary);
Objectives: 1 linear; 
Constraints:  176 linear;
Algebraic expressions:  1 min;
Logical expressions:  269 conditional (in)equalitie(s); 1540 not; 1345 or;

AMPL MP final model has 8029 variables (4515 integer, 3479 binary);
Objectives: 1 linear; 
Constraints:  3341 linear;


Running HiGHS 1.10.0 (git hash: fd86653): Copyright (c) 2025 HiGHS under MIT licence terms
MIP  has 3341 rows; 8029 cols; 27329 nonzeros; 3479 integer variables (3479 binary)
Coefficient ranges:
  Matrix [1e-01, 2e+04]
  Cost   [1e+00, 1e+00]
  Bound  [1e+00, 1e+04]
  RHS    [1e+00, 2e+04]
Presolving model
1800 rows, 1970 cols, 13372 nonzeros  0s
1524 rows, 1739 cols, 11368 nonzeros  0s

Solving MIP model with:
   1524 rows
   1739 cols (1708 binary, 0 integer, 0 implied int., 31 continuous)
   11368 nonzeros

Src: B => Branching; C => Central rounding; F => Feasibility pump; H => Heuristic; L => Sub-MIP;
     P => Empty MIP; R => Randomized rounding; S => Solve LP; T => Evaluate node; U => Unbounded;
     z => Trivial zero; l => Trivial lower; u => Trivial upper; p => Trivial point; X => User solution

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
Src  Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   6661.577577     -inf                 inf        0      0      0         0     0.5s
         0       0         0   0.00%   194.330916      -inf                 inf        0      0      7       472     0.5s
         0       0         0   0.00%   165.1028047     -inf                 inf     2021    106    215      3325     5.6s
 L       0       0         0   0.00%   164.2864891     27.98506821      487.05%     2307    120    259      4108     8.7s
 L       0       0         0   0.00%   164.0904875     71.37049922      129.91%     2312    122    259      8155     9.3s

Symmetry detection completed in 0.0s
Found 4 generator(s)

 L       0       0         0   0.00%   164.0904875     71.37050022      129.91%     2070     62    259      8629    11.5s
         6       0         1   3.12%   164.0904875     71.37050022      129.91%     2071     62    292     56623    17.9s
        13       1         4   7.42%   164.0904875     71.37050022      129.91%     2085     62    548    105765    23.3s
       209      31        72  14.38%   161.5512555     71.37050022      126.36%     2322    102   2432    142502    28.7s
       461      52       179  15.04%   161.5512555     71.37050022      126.36%     1848    100   5526    169493    33.8s
 L     611      68       245  15.05%   161.5512555     72.52622979      122.75%     1761     87   8453    178165    42.3s
       948      91       397  15.07%   161.5512555     72.52622979      122.75%     1736    103   8151    219237    47.3s

Restarting search from the root node
Model after restart has 1506 rows, 1721 cols (1690 bin., 0 int., 0 impl., 31 cont.), and 11141 nonzeros

      1215       0         0   0.00%   161.5512555     72.52622979      122.75%      103      0      0    244609    52.2s
      1215       0         0   0.00%   161.5512555     72.52622979      122.75%      103     47     16    245265    52.3s

Symmetry detection completed in 0.0s
Found 1 generator(s)

      1330      12        43   6.95%   161.5512555     72.52622979      122.75%     1166    114   1124    267088    57.4s
      1640      39       176   9.09%   161.5512555     72.52622979      122.75%     1403     95   4866    297142    62.4s
      1882      56       282  12.16%   161.5512555     72.52622979      122.75%     1437    104   8221    326262    67.5s
      2069      68       366  12.86%   161.5512555     72.52622979      122.75%     1379    123   9765    355010    72.6s
      2229      78       434  14.04%   161.5512555     72.52622979      122.75%     1945    139   9442    387334    77.7s
      2435      93       522  15.37%   161.5512555     72.52622979      122.75%     2162    128   8696    412998    82.8s

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
Src  Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

 T    2479      89       542  15.42%   161.5512555     89.72445048       80.05%     2175    128   9293    416682    83.5s
 T    2500      90       548  15.42%   161.5512555     89.8877421        79.73%     2179    128   9454    418176    83.8s
      2732     114       644  16.46%   161.5512555     89.8877421        79.73%     2409    172   9945    445187    88.8s
      2986     131       756  16.51%   161.5512555     89.8877421        79.73%     2413    125   9877    466773    93.9s
      3233     150       869  16.55%   161.5512555     89.8877421        79.73%     2829    121   8721    490341    98.9s
      3516     169       990  20.67%   161.5512555     89.8877421        79.73%     3499    139   9696    514476   104.0s
      3750     196      1086  23.33%   161.025181      89.8877421        79.14%     3765    114   7381    532437   109.0s
      3950     217      1169  23.94%   161.025181      89.8877421        79.14%     3883    136   8363    561214   114.1s
      4137     228      1258  23.97%   161.0239575     89.8877421        79.14%     3798    180   9467    583765   119.1s
      4322     240      1340  24.02%   161.0239575     89.8877421        79.14%     3848    203   7326    607279   124.1s
      4345     250      1349  24.03%   161.0239575     89.8877421        79.14%     3684    140   7602    644436   138.0s
      4542     263      1438  24.04%   161.0239575     89.8877421        79.14%     3714    145   9190    671231   143.0s
      4724     271      1511  24.06%   161.0239575     89.8877421        79.14%     2841    162   8176    700200   148.2s
      4785     290      1537  24.06%   161.0239575     89.8877421        79.14%     2769    141   9393    710564   150.0s

Solving report
  Status            Time limit reached
  Primal bound      89.8877420991
  Dual bound        161.023957479
  Gap               79.14% (tolerance: 0.01%)
  P-D integral      148.344174924
  Solution status   feasible
                    89.8877420991 (objective)
                    0 (bound viol.)
                    5.59659185448e-14 (int. viol.)
                    0 (row viol.)
  Timing            150.01 (total)
                    0.00 (presolve)
                    0.00 (solve)
                    0.00 (postsolve)
  Max sub-MIP depth 4
  Nodes             4785
  Repair LPs        0 (0 feasible; 0 iterations)
  LP iterations     710564 (total)
                    283463 (strong br.)
                    14059 (separation)
                    60177 (heuristics)
  Warning code 1 while solving with HiGHS
HiGHS 1.10.0: time limit, feasible solution; objective 89.8877421
710564 simplex iterations
4785 branching nodes
absmipgap=71.1362, relmipgap=0.79139

# Debug
def extract_cycles(arcs):
    visited = set()
    cycles = []

    # Create a mapping from i -> j
    arc_dict = {i: j for i, j in arcs}

    for start in arc_dict:
        if start not in visited:
            cycle = []
            current = start
            while current not in visited:
                visited.add(current)
                cycle.append(current)
                current = arc_dict.get(current)
                if current is None:
                    break  # Not a cycle or incomplete path
            # Close the cycle if it loops back to the start
            if cycle and current == cycle[0]:
                cycles.append(cycle)
    return cycles


# Debug


def print_loops(model):
    sol = model.var["x"].to_dict()
    V = model.set["V"]
    K = model.set["K"]
    for k in K:
        solk = [(i, j) for i in V for j in V if sol[i, j, k] >= 0.5]
        cycles = extract_cycles(solk)
        print(f"k = {k}:")
        for cycle in cycles:
            print(", ".join(map(str, cycle)) + f", {cycle[0]}")


m.display("MinProfit, TotalProfit")
print_loops(model=m)

MinProfit = 89.8877
TotalProfit = 671.369

k = 0:
0, 21, 2, 1, 14, 22, 20, 0
k = 1:
0, 15, 5, 3, 18, 7, 19, 0
k = 2:
0, 16, 6, 10, 24, 0
k = 3:
0, 23, 9, 17, 25, 0
k = 4:
0, 11, 13, 4, 8, 12, 0