# Cannot Get Optimal Solution with 16 nodes of VRP with Time Windows

I have been wondering if my code is wrong or something. But why can't I get the optimality with only 16 nodes (I have clustered it so the max now becomes 16 from a total of 54 nodes)? My model was based on VRP Time Windows by Paschos 2014. I've tried it by using Lingo and Gurobi. The formulation is this: $$\text { Minimize } \sum_{k=1}^{m} \sum_{\left(v_{i}, v_{j}\right) \in \mathcal{A}} c_{i j} x_{i j k}$$

Constraint: $$\begin{array}{ll} \sum_{k=1}^{m} \sum_{v_{j} \in \mathcal{V}} x_{i j k}=1 & \begin{array}{r} k=1\left(v_{i}, v_{j}\right) \in \mathcal{A} \\ v_{i} \in \mathcal{V} \backslash\left\{v_{0}\right\} \end{array} \\ \sum_{v_{i} \in \mathcal{V}} x_{i \ell k}=\sum_{v_{j} \in \mathcal{V}} x_{\ell j k} & v_{\ell} \in \mathcal{V} \backslash\left\{v_{0}\right\}, 1 \leqslant k \leqslant m \\ \sum_{v_{j} \in \mathcal{V} \backslash\left\{v_{0}\right\}} x_{0 j k}=1 & 1 \leqslant k \leqslant m \\ \sum_{v_{i} \in \mathcal{V} \backslash\left\{v_{0}\right\}} x_{i 0 k}=1 & 1 \leqslant k \leqslant m \\ \sum_{v_{i} \in \mathcal{V} \backslash\left\{v_{0}\right\}} \sum_{v_{j} \in \mathcal{V}} q_{i} x_{i j k} \leqslant Q & 1 \leqslant k \leqslant m \end{array}$$

\begin{aligned} &a_{i} \sum_{v_{j} \in \mathcal{V}} x_{i j k} \leqslant u_{i k} \leqslant b_{i} \sum_{v_{j} \in \mathcal{V}} x_{i j k}, \quad v_{i} \in \mathcal{V}, 1 \leqslant k \leqslant m \\ &u_{i k}+s_{i}+t_{i j}-u_{j k} \leqslant T\left(1-x_{i j k}\right), \quad\left(v_{i}, v_{j}\right) \in \mathcal{A}, 1 \leqslant k \leqslant m \\ &u_{i k} \geqslant 0, \quad v_{i} \in \mathcal{V}, 1 \leqslant k \leqslant m \end{aligned}

The weird thing is when I tried to solve it with 10 nodes, the lingo solver was able to return the optimal solution in under 1 minute. But when I tried with real data, even for 10 hours running, the solver is still running.

EDIT:

Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (win64)
Thread count: 2 physical cores, 4 logical processors, using up to 4 threads
Optimize a model with 168 rows, 768 columns and 2976 nonzeros
Model fingerprint: 0x460da4ce
Model has 630 general constraints
Variable types: 48 continuous, 720 integer (720 binary)
Coefficient statistics:
Matrix range     [1e+00, 8e+01]
Objective range  [1e+00, 5e+01]
Bounds range     [1e+00, 1e+00]
RHS range        [1e+00, 1e+03]
GenCon rhs range [1e+01, 4e+01]
GenCon coe range [1e+00, 1e+00]
Presolve added 1272 rows and 210 columns
Presolve time: 0.17s
Presolved: 1440 rows, 978 columns, 11136 nonzeros
Variable types: 258 continuous, 720 integer (720 binary)

Root relaxation: objective 2.700000e+01, 67 iterations, 0.03 seconds (0.00 work units)

Nodes    |    Current Node    |     Objective Bounds      |     Work
Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

0     0   27.00000    0   30          -   27.00000      -     -    0s
0     0   30.02381    0   37          -   30.02381      -     -    0s
0     0   35.00000    0   36          -   35.00000      -     -    0s
0     0   43.00000    0   31          -   43.00000      -     -    0s
H    0     0                     281.0000000   43.00000  84.7%     -    0s
0     0   43.00000    0   31  281.00000   43.00000  84.7%     -    0s
0     0   43.00000    0   31  281.00000   43.00000  84.7%     -    1s
0     0   43.00000    0   30  281.00000   43.00000  84.7%     -    1s
0     0   43.00000    0   30  281.00000   43.00000  84.7%     -    1s
0     0   43.00000    0   30  281.00000   43.00000  84.7%     -    1s
0     2   43.00000    0   26  281.00000   43.00000  84.7%     -    1s
H   29    35                     198.0000000   43.00000  78.3%  21.5    1s
H   88    99                     168.0000000   43.00000  74.4%  14.0    1s
*  119   123              50     127.0000000   43.00000  66.1%  12.4    1s
H  456   391                     125.0000000   43.00000  65.6%  10.3    2s
H  463   379                     119.0000000   43.00000  63.9%  10.4    2s
817   658   43.15451   22   27  119.00000   43.00000  63.9%  15.9    5s
H 1006   719                     118.0000000   43.00000  63.6%  16.5    5s
H 1038   692                     114.0000000   43.00000  62.3%  16.6    5s
H 1073   678                     113.0000000   43.00000  61.9%  16.6    6s
H 1705   899                     111.0000000   43.00000  61.3%  16.3    7s
3202  1880   84.97293   66   13  111.00000   43.48152  60.8%  15.5   10s
8845  6167   93.21239   72   10  111.00000   44.28152  60.1%  15.6   15s
14423  9889 infeasible   64       111.00000   45.12155  59.3%  14.7   20s
19007 13198   51.85990   38   12  111.00000   45.74524  58.8%  14.4   25s
20856 14247   67.91148   60   30  111.00000   46.00000  58.6%  14.5   31s
20886 14267   86.76100   43   54  111.00000   46.00000  58.6%  14.4   35s
20987 14336   46.00000   36   33  111.00000   46.00000  58.6%  14.8   40s
22171 14926  100.20962  115   23  111.00000   46.00000  58.6%  15.1   45s
26640 16889 infeasible   69       111.00000   47.53714  57.2%  15.0   50s
32818 19229   69.72771   59   26  111.00000   51.16970  53.9%  14.8   55s
38186 21010   59.97461   45   34  111.00000   53.16667  52.1%  14.9   60s
42489 22265   90.33333   49   14  111.00000   54.63747  50.8%  14.9   65s
47836 24148   89.49744   77   31  111.00000   56.00000  49.5%  14.9   70s
53826 25927   65.11477   60   34  111.00000   56.79788  48.8%  14.8   75s
59884 27894   78.77307  154   23  111.00000   57.41955  48.3%  14.9   80s
65719 30032  101.20672   60   24  111.00000   58.09387  47.7%  14.9   85s
71150 33633 infeasible  115       111.00000   58.72779  47.1%  14.9   90s
74887 36114   64.05556   52   27  111.00000   59.09387  46.8%  14.9   95s
79892 39307  109.30875   49   24  111.00000   59.51082  46.4%  14.9  100s
83720 41742   71.64684   51   22  111.00000   59.91729  46.0%  15.0  105s
88446 44919   98.25571   45   10  111.00000   60.10905  45.8%  15.0  110s
93592 48237   85.96847   48   24  111.00000   60.43274  45.6%  15.1  115s
98541 51522   64.60322   51   25  111.00000   60.72943  45.3%  15.1  120s
102845 54267   99.00000   65   16  111.00000   60.95160  45.1%  15.1  125s
107394 57673   69.52727   64   22  111.00000   61.00000  45.0%  15.2  130s
112444 61028   80.75000   78   18  111.00000   61.12951  44.9%  15.2  135s
117523 64288   67.74127   75   16  111.00000   61.26834  44.8%  15.3  140s
122741 67689   63.32275   65   19  111.00000   61.45337  44.6%  15.3  145s
127607 70950  101.00000   49   12  111.00000   61.57306  44.5%  15.4  150s
132756 74210   65.47593   42   31  111.00000   61.78555  44.3%  15.4  155s
136990 76979   92.93333   77   15  111.00000   61.93750  44.2%  15.4  160s
142582 80382  104.00000   76   17  111.00000   62.00000  44.1%  15.4  165s
147283 83875  109.09497   64   19  111.00000   62.00000  44.1%  15.5  170s
152753 87143  103.55460   60   10  111.00000   62.12951  44.0%  15.5  175s
158126 90549   93.94203   55   19  111.00000   62.23102  43.9%  15.5  180s
163038 93669     cutoff   60       111.00000   62.32400  43.9%  15.6  185s
166698 95951   77.50173   77   15  111.00000   62.40379  43.8%  15.6  190s
171727 98819  105.61448   56   22  111.00000   62.52864  43.7%  15.6  195s
176952 102393 infeasible   99       111.00000   62.64973  43.6%  15.6  200s
181955 105410   87.20043   45   11  111.00000   62.75701  43.5%  15.6  205s
187145 108539   87.80647   75   35  111.00000   62.89160  43.3%  15.7  210s
192246 111837 infeasible   79       111.00000   62.99708  43.2%  15.7  215s
197389 115062   70.75900   85   13  111.00000   63.00000  43.2%  15.7  220s
201353 117893   93.06245   75   31  111.00000   63.02439  43.2%  15.7  225s
206195 120822   63.47779   54   24  111.00000   63.09691  43.2%  15.7  230s
210899 123760  105.00546   61   14  111.00000   63.14964  43.1%  15.7  235s
216382 127092  105.00000   74   14  111.00000   63.22222  43.0%  15.7  240s
222038 130637   90.00000   68   12  111.00000   63.32452  43.0%  15.7  245s
226492 133185   90.91646   44   13  111.00000   63.37571  42.9%  15.7  250s
230755 136152   64.69033   64   20  111.00000   63.41868  42.9%  15.7  255s
235252 138714   63.86574   57   31  111.00000   63.49554  42.8%  15.8  260s
240283 141816 infeasible   50       111.00000   63.55502  42.7%  15.7  265s
246089 145454   79.49498   48   18  111.00000   63.64156  42.7%  15.7  270s
250670 148171   79.08612   52   12  111.00000   63.70720  42.6%  15.7  275s
255030 150974   80.73087   37   15  111.00000   63.75478  42.6%  15.8  280s
258930 153198   94.81954   45   16  111.00000   63.79969  42.5%  15.8  285s
263816 156010   76.15049   45   25  111.00000   63.88839  42.4%  15.8  290s
268782 158824   93.50000   49    8  111.00000   63.95633  42.4%  15.8  295s
271310 160524   84.00000   52   15  111.00000   64.00000  42.3%  15.8  300s
274848 162583   66.86667   34   27  111.00000   64.00000  42.3%  15.8  305s
278496 165062 infeasible   87       111.00000   64.00000  42.3%  15.8  310s
281970 167381     cutoff   59       111.00000   64.01453  42.3%  15.8  315s
285890 169888   97.42763   61   20  111.00000   64.06042  42.3%  15.9  320s

• Can you show the logs of the solvers? Commented Jul 21, 2022 at 10:14
• sure, i will add the logs in edit Commented Jul 21, 2022 at 11:14
• Since the vehicles have the same capacities, there is no need to use a three-index formulation (which are known to be difficult to solve due to e.g. symmetry). Hence, try to formulate your problem using $x_{ij}$-variables instead of $x_{ijk} & - variables. – Sune Commented Jul 21, 2022 at 12:04 • homonegenous fleet, i see. Can you refer me the two indices time windows? Commented Jul 21, 2022 at 12:14 • Let$u_i$be the time a vehicle starts service at node$i$. Then$a_i\leq u_i\leq b_i$and$u_i - u_j +s_i +t_{ij} \leq M(1-x_{ij})$for all customers$i$and$j\$. Make sure to handle the "first customer after depot" case as well.
– Sune
Commented Jul 21, 2022 at 13:11

It is well-known that MIP solvers are relatively inefficient for solving VRPs. If you really want to obtain (and prove) an optimal solution, you should use branch-cut-and-price solvers like VRPSolver (https://vrpsolver.math.u-bordeaux.fr). If you prefer C++, there is BaPCod with VRPSolver extension (https://bapcod.math.u-bordeaux.fr/), it has the VRPTW demo which you can use. There is also GCG (https://gcg.or.rwth-aachen.de/), which is less efficient for VRPs, but still solving an instance with 16 nodes will be very fast with it. If it is sufficient to get just a good solution, there are heuristic solvers, like LocalSolver, OR-Tools, Optaplanner, as suggested by others.

• how to implement this? i have a basic python and zero knowledge julia Commented Jul 25, 2022 at 9:38
• Sorry, only Julia or C++ (at the moment). At the same time, there is the VRPTW demo already available, so no coding is required (if the datafile format is the same). Commented Jul 26, 2022 at 11:08
• There is also VRPY package in Python (github.com/Kuifje02/vrpy). It does not guarantee to return an optimal solution, but in some cases it will. May be it will do for such a small instance. Commented Jul 26, 2022 at 11:11

If you are sure that you don't have any bugs or flows in your code, I would recommend more strong solvers like CPLEX, Gurobi, Localsolver (those have free academic licenses) or at the very least Google OR Tools with its free license SCIP solver or its specialized Constraint Programming Solver which has a special VRPTW optimizer.

• i did try with gurobi, but the result is the same. The solver find difficult to get the optimality. Does the time windows width become more complex as they are getting bigger? I have time windows of 0 for all the earliest and 480 (in minutes) for the lastest. Commented Jul 21, 2022 at 10:55
• I'm not dive down into localsolver, but when i try using solomon instance, the result is not very good and it still difficult to find the optimality Commented Jul 21, 2022 at 10:56
• I think i'm down to try other solver, because the deadline is so much closer now Commented Jul 21, 2022 at 10:58
• Hi @OneAttack. I am surprised by your feedback. Note that LocalSolver is not a MILP solver. The modeling API is different, between mixed-integer linear programming and constraint programming. The CVRPTW model is given here: localsolver.com/docs/last/exampletour/vrptw.html. I would be surprised that you don't get near-optimal solutions in seconds on such small instances. You can check a benchmark here: localsolver.com/benchmark/…. As said above, LocalSolver is free for academics. Commented Jul 21, 2022 at 19:13

And yet another open source option (python): vrpy.

An example with time windows (borrowed from the excellent ortools module) is presented in the docs here:

from networkx import DiGraph, from_numpy_matrix, relabel_nodes, set_node_attributes
from numpy import array

# Distance matrix
DISTANCES = [
[0,548,776,696,582,274,502,194,308,194,536,502,388,354,468,776,662,0], # from Source
[0,0,684,308,194,502,730,354,696,742,1084,594,480,674,1016,868,1210,548],
[0,684,0,992,878,502,274,810,468,742,400,1278,1164,1130,788,1552,754,776],
[0,308,992,0,114,650,878,502,844,890,1232,514,628,822,1164,560,1358,696],
[0,194,878,114,0,536,764,388,730,776,1118,400,514,708,1050,674,1244,582],
[0,502,502,650,536,0,228,308,194,240,582,776,662,628,514,1050,708,274],
[0,730,274,878,764,228,0,536,194,468,354,1004,890,856,514,1278,480,502],
[0,354,810,502,388,308,536,0,342,388,730,468,354,320,662,742,856,194],
[0,696,468,844,730,194,194,342,0,274,388,810,696,662,320,1084,514,308],
[0,742,742,890,776,240,468,388,274,0,342,536,422,388,274,810,468,194],
[0,1084,400,1232,1118,582,354,730,388,342,0,878,764,730,388,1152,354,536],
[0,594,1278,514,400,776,1004,468,810,536,878,0,114,308,650,274,844,502],
[0,480,1164,628,514,662,890,354,696,422,764,114,0,194,536,388,730,388],
[0,674,1130,822,708,628,856,320,662,388,730,308,194,0,342,422,536,354],
[0,1016,788,1164,1050,514,514,662,320,274,388,650,536,342,0,764,194,468],
[0,868,1552,560,674,1050,1278,742,1084,810,1152,274,388,422,764,0,798,776],
[0,1210,754,1358,1244,708,480,856,514,468,354,844,730,536,194,798,0,662],
[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], # from Sink
]

TRAVEL_TIMES = [
[0, 6, 9, 8, 7, 3, 6, 2, 3, 2, 6, 6, 4, 4, 5, 9, 7, 0],  # from source
[0, 0, 8, 3, 2, 6, 8, 4, 8, 8, 13, 7, 5, 8, 12, 10, 14, 6],
[0, 8, 0, 11, 10, 6, 3, 9, 5, 8, 4, 15, 14, 13, 9, 18, 9, 9],
[0, 3, 11, 0, 1, 7, 10, 6, 10, 10, 14, 6, 7, 9, 14, 6, 16, 8],
[0, 2, 10, 1, 0, 6, 9, 4, 8, 9, 13, 4, 6, 8, 12, 8, 14, 7],
[0, 6, 6, 7, 6, 0, 2, 3, 2, 2, 7, 9, 7, 7, 6, 12, 8, 3],
[0, 8, 3, 10, 9, 2, 0, 6, 2, 5, 4, 12, 10, 10, 6, 15, 5, 6],
[0, 4, 9, 6, 4, 3, 6, 0, 4, 4, 8, 5, 4, 3, 7, 8, 10, 2],
[0, 8, 5, 10, 8, 2, 2, 4, 0, 3, 4, 9, 8, 7, 3, 13, 6, 3],
[0, 8, 8, 10, 9, 2, 5, 4, 3, 0, 4, 6, 5, 4, 3, 9, 5, 2],
[0, 13, 4, 14, 13, 7, 4, 8, 4, 4, 0, 10, 9, 8, 4, 13, 4, 6],
[0, 7, 15, 6, 4, 9, 12, 5, 9, 6, 10, 0, 1, 3, 7, 3, 10, 6],
[0, 5, 14, 7, 6, 7, 10, 4, 8, 5, 9, 1, 0, 2, 6, 4, 8, 4],
[0, 8, 13, 9, 8, 7, 10, 3, 7, 4, 8, 3, 2, 0, 4, 5, 6, 4],
[0, 12, 9, 14, 12, 6, 6, 7, 3, 3, 4, 7, 6, 4, 0, 9, 2, 5],
[0, 10, 18, 6, 8, 12, 15, 8, 13, 9, 13, 3, 4, 5, 9, 0, 9, 9],
[0, 14, 9, 16, 14, 8, 5, 10, 6, 5, 4, 10, 8, 6, 2, 9, 0, 7],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # from sink
]

# Time windows (key: node, value: lower/upper bound)
TIME_WINDOWS_LOWER = {0: 0, 1: 7, 2: 10, 3: 16, 4: 10, 5: 0, 6: 5, 7: 0, 8: 5, 9: 0, 10: 10, 11: 10, 12: 0, 13: 5, 14: 7, 15: 10, 16: 11,}
TIME_WINDOWS_UPPER = {1: 12, 2: 15, 3: 18, 4: 13, 5: 5, 6: 10, 7: 4, 8: 10, 9: 3, 10: 16, 11: 15, 12: 5, 13: 10, 14: 8, 15: 15, 16: 15,}

# Transform distance matrix into DiGraph
A = array(DISTANCES, dtype=[("cost", int)])
G_d = from_numpy_matrix(A, create_using=DiGraph())

# Transform time matrix into DiGraph
A = array(TRAVEL_TIMES, dtype=[("time", int)])
G_t = from_numpy_matrix(A, create_using=DiGraph())

# Merge
G = compose(G_d, G_t)

# Set time windows
set_node_attributes(G, values=TIME_WINDOWS_LOWER, name="lower")
set_node_attributes(G, values=TIME_WINDOWS_UPPER, name="upper")

# The VRP is defined and solved
prob = VehicleRoutingProblem(G, time_windows=True)
prob.solve()

• 'I have been in the OR field for 10 years. Still lovin it.' --> what happened to questions to all your answers?
– BCLC
Commented Jul 22, 2022 at 16:06
• aha you got me ;) Commented Jul 24, 2022 at 18:51

As the comment above explains, the Capacitated Vehicle Routing Problem with Time Windows (CVRPTW) can be modeled compactly by following a list-based modeling approach instead of the classical Boolean modeling approach, as you presented in your question. This modeling approach offered by Hexaly is different from traditional MILP solvers. If you wish to get good results using Hexaly on problems like vehicle routing, job shop scheduling, or bin packing, please use this set or list-based modeling approach.

Note that LocalSolver is commercial software. Nevertheless, it is free for faculty and students.

Below is the code snippet to model the CVRPTW by using Hexaly Modeler:

function model() {
customersSequences[k in 1..nbTrucks] <- list(nbCustomers);

// All customers must be visited by the trucks
constraint partition[k in 1..nbTrucks](customersSequences[k]);

for[k in 1..nbTrucks] {
local sequence <- customersSequences[k];
local c <- count(sequence);

// A truck is used if it visits at least one customer
truckUsed[k] <- c > 0;

// The quantity needed in each route must not exceed the truck capacity
routeQuantity[k] <- sum(0..c-1, i => demands[sequence[i]]);
constraint routeQuantity[k] <= truckCapacity;

endTime[k] <- array(0..c-1, (i, prev) => max(earliestStart[sequence[i]],
i == 0 ? distanceWarehouse[sequence[0]] :
prev + distanceMatrix[sequence[i-1]][sequence[i]])
+ serviceTime[sequence[i]]);

homeLateness[k] <- truckUsed[k] ?
max(0, endTime[k][c - 1] + distanceWarehouse[sequence[c - 1]] - maxHorizon) :
0;

// Distance traveled by truck k
routeDistances[k] <- sum(1..c-1,
i => distanceMatrix[sequence[i-1]][sequence[i]]) + (truckUsed[k] ?
(distanceWarehouse[sequence[0]] + distanceWarehouse[sequence[c-1]]) :
0);

lateness[k] <- homeLateness[k] + sum(0..c-1,
i => max(0, endTime[k][i] - latestEnd[sequence[i]]));
}

// Total lateness, must be 0 for a solution to be valid
totalLateness <- sum[k in 1..nbTrucks](lateness[k]);

nbTrucksUsed <- sum[k in 1..nbTrucks](truckUsed[k]);

// Total distance traveled (convention in Solomon's instances is to round to 2 decimals)
totalDistance <- round(100 * sum[k in 1..nbTrucks](routeDistances[k])) / 100;

minimize totalLateness;
minimize nbTrucksUsed;
minimize totalDistance;
}


For the complete ready-to-use code in Hexaly Modeler, Python, Java, C#, or C++, please have a look at https://www.hexaly.com/docs/last/exampletour/vrptw.html

Another alternative for VRPTW is OptaPlanner (Java) or OptaPy (Python). It's open source (Apache license, so 100% free) and used across the globe extensively, including in Fortune 500 companies that use it for VRPTW cases of 10 000+ vehicles.

Example source code to assign each customer to a vehicle:

Model

@PlanningEntity
public class Vehicle {

private long id;
private int capacity;
private Depot depot;

@PlanningListVariable(valueRangeProviderRefs = "customerRange")
private List<Customer> customerList;

...
}

public class Customer {

private long id;
private Location location;
private int demand;

...
}


Constraints

public class VehicleRoutingConstraintProvider implements ConstraintProvider {

@Override
public Constraint[] defineConstraints(ConstraintFactory factory) {
return new Constraint[] {
vehicleCapacity(factory),
totalDistance(factory),
};
}

// ************************************************************************
// Hard constraints
// ************************************************************************

protected Constraint vehicleCapacity(ConstraintFactory factory) {
return factory.forEach(Vehicle.class)
.filter(vehicle -> vehicle.getTotalDemand() > vehicle.getCapacity())
.penalizeLong(
"vehicleCapacity",
HardSoftLongScore.ONE_HARD,
vehicle -> vehicle.getTotalDemand() - vehicle.getCapacity());
}

// ************************************************************************
// Soft constraints
// ************************************************************************

protected Constraint totalDistance(ConstraintFactory factory) {
return factory.forEach(Vehicle.class)
.penalizeLong(
"distanceFromPreviousStandstill",
HardSoftLongScore.ONE_SOFT,
Vehicle::getTotalDistanceMeters);
}
}


Generic MILP solvers cannot solve VRPTWs to optimality. You need to use methods specifically designed for VRPTW. Column generation methods (branch-cut-and-price) are the most effective for this class of problems.

• Do you know a good starting point how to use column generation? Commented Jul 26, 2022 at 12:14
• This is a great tutorial on the column generation approach to solve VRPs: link.springer.com/article/10.1007/s10288-010-0130-z Commented Jul 26, 2022 at 14:11