• Docs »
  • Tutorials »
  • Modeling Routing Problems with List Decisions

Modeling Routing Problems with List Decisions

List decision variables are a very powerful modeling feature for various problems where collections of items have to be ordered in an optimized fashion. For instance, routing, scheduling, production planning, or even assignment problems are typical use cases. Indeed, list-based models for such problems are very compact and straightforward. Not only are they easier to write than classic MILP models, they also allow Hexaly to scale well and achieve state-of-the-art results even on very large instances. For more information, check out our routing benchmarks page.

Starting with the Traveling Salesman Problem, this page illustrates the good modeling practices for routing problems through guided step-by-step examples.

The Traveling Salesman Problem (TSP)

The Traveling Salesman Problem (TSP) is defined as follows. Given a set of n cities and distances between each pair of cities, find a roundtrip of minimal total length visiting each city exactly once. The distance from city i to city j and the distance from city j to city i may be different.

While it is possible to write a linear program for the Traveling Salesman Problem (TSP), classic MILP models for this problem are not so straightforward. Indeed, enforcing connectivity is difficult, and classic MILP formulations rely on an exponential number of subtour-elimination constraints or auxiliary variables to ensure that the solution forms a single tour, leading to models that are substantially more intricate than the problem statement itself. Moreover, such models typically require a large number of binary decision variables, making it harder for solvers to efficiently solve large instances.

On the contrary, the natural modeling of the Traveling Salesman Problem (TSP) only suggests finding an optimal permutation of the cities, which is precisely what list decision variables are for. Indeed, a list decision of domain size n represents a subpermutation of {0, ..., n-1}.

The Hexaly model for the Traveling Salesman Problem (TSP) starts with the decisions definition. It uses a single list decision, named cities, representing the permutation of cities: cities[i] represents the index of the city in position i in the tour:

cities <- list(n);

cities = model.list(nb_cities)

cities = model.listVar(nbCities);

cities = model.List(nbCities);

cities = model.listVar(nbCities);

Then, we need to add a constraint ensuring that all cities are visited:

constraint count(cities) == n;

model.constraint(model.count(cities) == nb_cities)

model.constraint(model.count(cities) == nbCities);

model.Constraint(model.Count(cities) == nbCities);

model.constraint(model.eq(model.count(cities), nbCities));

The count operator returns the number of elements in the list (size of the tour). Since the list decision represents a permutation, every element inside cities is unique: ensuring that it contains n elements ensures that every city is visited exactly once. This single constraint guarantees that the solution is valid and forms a single connected tour of all cities.

The only remaining step consists in writing the objective function. We compute the total traveled distance, which is the sum of pairwise distances from one city to the next, and from the last city back to the first one. In the following code, totalDistance is not a decision variable but an intermediate expression: its value is directly computed from the value of the list decision and the input data. Then, we indicate that this expression is an objective function to minimize:

totalDistance <- sum(1...n, i => distance[cities[i-1]][cities[i]])
        + distance[cities[n-1]][cities[0]];
minimize totalDistance;

dist_lambda = model.lambda_function(lambda i:
                                    model.at(dist_matrix, cities[i - 1], cities[i]))
total_distance = model.sum(model.range(1, nb_cities), dist_lambda) \
    + model.at(dist_matrix, cities[nb_cities - 1], cities[0])
model.minimize(total_distance)

HxExpression distLambda =
    model.createLambdaFunction([&](HxExpression i) { return model.at(distMatrix, cities[i - 1], cities[i]); });
totalDistance = model.sum(model.range(1, nbCities), distLambda) + model.at(distMatrix, cities[nbCities - 1], cities[0]);
model.minimize(totalDistance);

HxExpression distLambda = model.LambdaFunction(i => distMatrix[cities[i - 1], cities[i]]);
totalDistance = model.Sum(model.Range(1, nbCities), distLambda) + distMatrix[cities[nbCities - 1], cities[0]];
model.Minimize(totalDistance);

HxExpression distLambda = model
    .lambdaFunction(i -> model.at(distMatrix, model.at(cities, model.sub(i, 1)), model.at(cities, i)));
totalDistance = model.sum(model.sum(model.range(1, nbCities), distLambda),
    model.at(distMatrix, model.at(cities, nbCities - 1), model.at(cities, 0)));
model.minimize(totalDistance);

This formulation uses a variadic sum operator and a lambda expression. The first line is equivalent to the following mathematical notation:

\[\sum_{i = 0}^{n-1}{d_{c_{i-1}, c_i}}\]

with distance[cities[i-1]][cities[i]] the distance between the cities in positions i-1 and i in the tour according to the value of the list decision variable.

The complete Hexaly model for the Traveling Salesman Problem (TSP) is then:

function model() {
    cities <- list(n);
    constraint count(cities) == n;
    totalDistance <- sum(1...n, i => distance[cities[i-1]][cities[i]])
            + distance[cities[n-1]][cities[0]];
    minimize totalDistance;
}

cities = model.list(nb_cities)

model.constraint(model.count(cities) == nb_cities)

dist_matrix = model.array(dist_matrix_data)
dist_lambda = model.lambda_function(lambda i:
                                    model.at(dist_matrix, cities[i - 1], cities[i]))
total_distance = model.sum(model.range(1, nb_cities), dist_lambda) \
    + model.at(dist_matrix, cities[nb_cities - 1], cities[0])
model.minimize(total_distance)

cities = model.listVar(nbCities);

model.constraint(model.count(cities) == nbCities);

HxExpression distMatrix = model.array();
for (int i = 0; i < nbCities; ++i) {
    HxExpression row = model.array(distMatrixData[i].begin(), distMatrixData[i].end());
    distMatrix.addOperand(row);
}
HxExpression distLambda =
    model.createLambdaFunction([&](HxExpression i) { return model.at(distMatrix, cities[i - 1], cities[i]); });
totalDistance = model.sum(model.range(1, nbCities), distLambda) + model.at(distMatrix, cities[nbCities - 1], cities[0]);
model.minimize(totalDistance);

cities = model.List(nbCities);

model.Constraint(model.Count(cities) == nbCities);

HxExpression distMatrix = model.Array(distMatrixData);
HxExpression distLambda = model.LambdaFunction(i => distMatrix[cities[i - 1], cities[i]]);
totalDistance = model.Sum(model.Range(1, nbCities), distLambda) + distMatrix[cities[nbCities - 1], cities[0]];
model.Minimize(totalDistance);

cities = model.listVar(nbCities);

model.constraint(model.eq(model.count(cities), nbCities));

HxExpression distMatrix = model.array(distMatrixData);
HxExpression distLambda = model
    .lambdaFunction(i -> model.at(distMatrix, model.at(cities, model.sub(i, 1)), model.at(cities, i)));
totalDistance = model.sum(model.sum(model.range(1, nbCities), distLambda),
    model.at(distMatrix, model.at(cities, nbCities - 1), model.at(cities, 0)));
model.minimize(totalDistance);

The whole code for solving the Traveling Salesman Problem (TSP) with Hexaly, along with a few instances, is also available as a code template, with an associated TSP benchmark.

The Capacitated Vehicle Routing Problem (CVRP)

The Capacitated Vehicle Routing Problem (CVRP) is an extension of the Traveling Salesman Problem (TSP) with several vehicles instead of one.

In the Capacitated Vehicle Routing Problem (CVRP), a fleet of delivery vehicles with uniform capacity must service customers with known demand for a single commodity. The vehicles start and end their routes at a common depot. Each customer must be served by exactly one vehicle, and the total demand served by each vehicle must not exceed its capacity. The objectives are to minimize the fleet size and the total traveled distance.

Similarly to the Traveling Salesman Problem (TSP), the Capacitated Vehicle Routing Problem (CVRP) uses list decision variables representing the order of customers visited by each truck:

customersSequences[k in 0...nbTrucks] <- list(nbCustomers);

customers_sequences = [model.list(nb_customers) for _ in range(nb_trucks)]

for (int k = 0; k < nbTrucks; ++k) {
    customersSequences[k] = model.listVar(nbCustomers);
}

for (int k = 0; k < nbTrucks; ++k)
    customersSequences[k] = model.List(nbCustomers);

for (int k = 0; k < nbTrucks; ++k)
    customersSequences[k] = model.listVar(nbCustomers);

For each truck k, customersSequences[k] is a list decision variable of domain size nbCustomers representing the customers served by truck k in visiting order: customersSequences[k][i] is the index of the customer served in position i in truck k’s tour.

All customers need to be served exactly once, by either one of the trucks. We use the partition operator to enforce this constraint:

constraint partition[k in 0...nbTrucks](customersSequences[k]);

model.constraint(model.partition(customers_sequences))

model.constraint(model.partition(customersSequences.begin(), customersSequences.end()));

model.Constraint(model.Partition(customersSequences));

model.constraint(model.partition(customersSequences));

The partition operator, when applied to a series of lists of domain size n, returns 1 if every element between 0 and n-1 is present in exactly one of the lists. Therefore, this constraint ensures that the index of each customer is present in exactly one of the lists, meaning that each customer is served by exactly one of the trucks.

To enforce the capacity constraints on the trucks, we need to compute the total demand of customers served by each truck. For that, we use the variadic form of the sum operator and a lambda function returning the demand of a customer:

local sequence <- customersSequences[k];
routeQuantity <- sum(sequence, j => demands[j]);
constraint routeQuantity <= truckCapacity;

sequence = customers_sequences[k]
demand_lambda = model.lambda_function(lambda j: demands[j])
route_quantity = model.sum(sequence, demand_lambda)
model.constraint(route_quantity <= truck_capacity)

HxExpression sequence = customersSequences[k];
HxExpression demandLambda =
        model.createLambdaFunction([&](HxExpression j) { return demands[j]; });
HxExpression routeQuantity = model.sum(sequence, demandLambda);
model.constraint(routeQuantity <= truckCapacity);

HxExpression sequence = customersSequences[k];
HxExpression demandLambda = model.LambdaFunction(j => demands[j]);
HxExpression routeQuantity = model.Sum(sequence, demandLambda);
model.Constraint(routeQuantity <= truckCapacity);

HxExpression sequence = customersSequences[k];
HxExpression demandLambda = model.lambdaFunction(j -> model.at(demands, j));
HxExpression routeQuantity = model.sum(sequence, demandLambda);
model.constraint(model.leq(routeQuantity, truckCapacity));

This expression loops over all the elements in customersSequences[k] to sum the associated demands. Note that the number of terms in this sum varies during the search (it is equal to the number of customers served by truck k). This feature of Hexaly Optimizer’s modeling formalism is particularly powerful, since being able to apply n-ary operators to a variable number of terms can make the models much more compact, and lead to better performance overall on large instances.

The computation of the distance traveled by each truck is similar to the objective function in the Traveling Salesman Problem:

local c <- count(sequence);
routeDistances[k] <- sum(1...c, i => distanceMatrix[sequence[i - 1]][sequence[i]])
            + (c > 0 ? (distanceDepot[sequence[0]] + distanceDepot[sequence[c - 1]]) : 0);

c = model.count(sequence)
dist_lambda = model.lambda_function(lambda i:
                                    model.at(dist_matrix, sequence[i - 1], sequence[i]))
dist_routes[k] = model.sum(model.range(1, c), dist_lambda) \
    + model.iif(c > 0,
                dist_depot[sequence[0]] + dist_depot[sequence[c - 1]],
                0)

HxExpression c = model.count(sequence);
HxExpression distLambda = model.createLambdaFunction(
    [&](HxExpression i) { return model.at(distMatrix, sequence[i - 1], sequence[i]); });
distRoutes[k] = model.sum(model.range(1, c), distLambda) +
                model.iif(c > 0, distDepot[sequence[0]] + distDepot[sequence[c - 1]], 0);

HxExpression c = model.Count(sequence);
HxExpression distLambda = model.LambdaFunction(
    i => distMatrix[sequence[i - 1], sequence[i]]
);
distRoutes[k] =
    model.Sum(model.Range(1, c), distLambda)
    + model.If(c > 0, distDepot[sequence[0]] + distDepot[sequence[c - 1]], 0);

HxExpression c = model.count(sequence);
HxExpression distLambda = model
    .lambdaFunction(i -> model.at(distMatrix, model.at(sequence, model.sub(i, 1)), model.at(sequence, i)));
distRoutes[k] = model.sum(model.sum(model.range(1, c), distLambda),
    model.iif(model.gt(c, 0), model.sum(model.at(distDepot, model.at(sequence, 0)),
        model.at(distDepot, model.at(sequence, model.sub(c, 1)))), 0));

If the truck does not serve any customer, the traveled distance is 0. Else, it is equal to the sum of distances from each customer to the next one (there again, note that the number of terms in the sum depends on the size of the list decision associated with the truck, and therefore varies during the search) plus the distance from the depot to the first customer and the distance from the last customer back to the depot.

Finally, we declare two objective functions, which are optimized in lexicographic order: objectives declared earlier take priority over those declared later. Here, we first minimize the number of trucks used, then minimize the total traveled distance. A solution using fewer trucks is therefore preferred even if its total traveled distance is greater. Among solutions using the same number of trucks, the one with the shortest total traveled distance is preferred.

nbTrucksUsed <- sum[k in 0...nbTrucks](trucksUsed[k]);
totalDistance <- sum[k in 0...nbTrucks](routeDistances[k]);

minimize nbTrucksUsed;
minimize totalDistance;

nb_trucks_used = model.sum(trucks_used)
total_distance = model.sum(dist_routes)

model.minimize(nb_trucks_used)
model.minimize(total_distance)

nbTrucksUsed = model.sum(trucksUsed.begin(), trucksUsed.end());
totalDistance = model.sum(distRoutes.begin(), distRoutes.end());

model.minimize(nbTrucksUsed);
model.minimize(totalDistance);

nbTrucksUsed = model.Sum(trucksUsed);
totalDistance = model.Sum(distRoutes);

model.Minimize(nbTrucksUsed);
model.Minimize(totalDistance);

nbTrucksUsed = model.sum(trucksUsed);
totalDistance = model.sum(distRoutes);

model.minimize(nbTrucksUsed);
model.minimize(totalDistance);

The complete Hexaly model for the Capacitated Vehicle Routing Problem (CVRP), along with a few instances, is available as a code template, with an associated CVRP benchmark.

The Capacitated Vehicle Routing Problem with Time Windows (CVRPTW)

The Capacitated Vehicle Routing Problem with Time Windows (CVRPTW) is an extension of the Capacitated Vehicle Routing Problem (CVRP) in which the customers must be served within opening hours (time windows).

These time windows can be difficult to respect in real-life applications. In practice, rather than concluding that no solution exists, it is often preferable to return a solution violating the time windows as little as possible. Consequently, the time windows are handled as a first-priority objective rather than hard constraints. In case of early arrival, the truck has to wait for the customer’s opening hour. In case of late arrival, we measure and penalize the lateness. The first objective of the problem then consists in minimizing the total lateness among all customers. A solution is considered feasible when this cumulated lateness is zero.

To compute the lateness, we first need to compute the end times of each truck’s visits. For each truck and each customer, the arrival time is the maximum of:

  • the end time of the previous visit plus the travel time (equivalent to the distance in this example). For the first visit, it is the travel time from the depot (assuming the truck leaves the depot at time 0)
  • the earliest allowed arrival time for this customer.

The end time is merely this arrival time plus the service time for this customer. The arrival to the depot at the end of the tour is the end time of the last visit plus the traveling time from this point to the depot.

The end times of each truck are computed using the array operator and a lambda function to create a recursive array: we allow a recursive definition of elements of this array by using a second argument in the function, containing the evaluation of the function on the previous element of the range:

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

end_time_lambda = model.lambda_function(
    lambda i, prev:
        model.max(
            earliest[sequence[i]],
            model.iif(
                i == 0, dist_depot[sequence[0]],
                prev + model.at(dist_matrix, sequence[i - 1], sequence[i])))
        + service_time[sequence[i]])

end_time[k] = model.array(model.range(0, c), end_time_lambda, 0)

HxExpression endTimeLambda = model.createLambdaFunction([&](HxExpression i, HxExpression prev) {
    return model.max(earliest[sequence[i]],
                     model.iif(i == 0, distDepot[sequence[0]],
                               prev + model.at(distMatrix, sequence[i - 1], sequence[i]))) +
           serviceTime[sequence[i]];
});
endTime[k] = model.array(model.range(0, c), endTimeLambda, 0);

HxExpression endTimeLambda = model.LambdaFunction(
    (i, prev) =>
        model.Max(
            earliest[sequence[i]],
            model.If(
                i == 0,
                distDepot[sequence[0]],
                prev + distMatrix[sequence[i - 1], sequence[i]]
            )
        ) + serviceTime[sequence[i]]
);
endTime[k] = model.Array(model.Range(0, c), endTimeLambda, 0);

HxExpression endTimeLambda = model.lambdaFunction((i, prev) -> model.sum(
    model.max(model.at(earliest, model.at(sequence, i)),
        model.sum(model.iif(model.eq(i, 0), model.at(distDepot, model.at(sequence, 0)),
            model.sum(prev, model.at(distMatrix, model.at(sequence, model.sub(i, 1)), model.at(sequence, i)))))),
    model.at(serviceTime, model.at(sequence, i))));

endTime[k] = model.array(model.range(0, c), endTimeLambda, 0);

Here, prev corresponds to the previous element in the array: the end time of the previous visit. From the end times, we can then easily compute the lateness associated with each truck k:

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

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

home_lateness[k] = model.iif(
    trucks_used[k],
    model.max(
        0,
        end_time[k][c - 1] + dist_depot[sequence[c - 1]] - max_horizon),
    0)

late_lambda = model.lambda_function(
    lambda i: model.max(0, end_time[k][i] - latest[sequence[i]]))
lateness[k] = home_lateness[k] + model.sum(model.range(0, c), late_lambda)

homeLateness[k] =
    model.iif(trucksUsed[k], model.max(0, endTime[k][c - 1] + distDepot[sequence[c - 1]] - maxHorizon), 0);

HxExpression lateLambda = model.createLambdaFunction(
    [&](HxExpression i) { return model.max(0, endTime[k][i] - latest[sequence[i]]); });
lateness[k] = homeLateness[k] + model.sum(model.range(0, c), lateLambda);

homeLateness[k] = model.If(
    trucksUsed[k],
    model.Max(0, endTime[k][c - 1] + distDepot[sequence[c - 1]] - maxHorizon),
    0
);

HxExpression lateLambda = model.LambdaFunction(
    i => model.Max(endTime[k][i] - latest[sequence[i]], 0)
);
lateness[k] = homeLateness[k] + model.Sum(model.Range(0, c), lateLambda);

HxExpression theEnd = endTime[k];

homeLateness[k] = model.iif(trucksUsed[k],
    model.max(0,
        model.sum(model.at(theEnd, model.sub(c, 1)), model.sub(model.at(distDepot, model.at(sequence, model.sub(c, 1))), maxHorizon))),
    0);

HxExpression lateLambda = m
    .lambdaFunction(i -> model.max(model.sub(model.at(theEnd, i), model.at(latest, model.at(sequence, i))), 0));
lateness[k] = model.sum(homeLateness[k], model.sum(model.range(0, c), lateLambda));

We can then minimize the cumulated lateness:

totalLateness <- sum[k in 0...nbTrucks](lateness[k]);
minimize totalLateness;

total_lateness = model.sum(lateness)
model.minimize(total_lateness)

totalLateness = model.sum(lateness.begin(), lateness.end());
model.minimize(totalLateness);

totalLateness = model.Sum(lateness);
model.Minimize(totalLateness);

totalLateness = model.sum(lateness);
model.minimize(totalLateness);

The complete Hexaly model for the Capacitated Vehicle Routing Problem with Time Windows (CVRPTW), along with a few instances, is available as a code template, with an associated CVRPTW benchmark.

Further reading: discover our code templates on routing problems

Pickup and Delivery constraints

When considering a routing problem with pickup and delivery constraints, the customers are paired: in each pair, one customer corresponds to a pickup point (with positive demand), and the other corresponds to a delivery point (with negative demand). Paired customers must be served by the same truck, and each pickup must be treated before its associated delivery.

The complete Hexaly model for the Pickup and Delivery Problem with Time Windows (PDPTW), along with a few instances, is available as a code template, with an associated PDPTW benchmark.

Prize-Collecting Routing Problems

In prize-collecting variants of routing problems, visiting all customers is unnecessary, but there is a minimum demand to satisfy. Furthermore, each visited customer provides a certain prize, and the total collected prize must be maximized.

The complete Hexaly model for the Prize-Collecting Vehicle Routing Problem (PCVRP), along with a few instances, is available as a code template.

Split Delivery

In split delivery routing problems, the goods delivered to the customers can be split among several trucks. In other words, each customer can be served by multiple vehicles and the total quantity served to each customer must be superior or equal to their demand.

The complete Hexaly model for the Split Delivery Vehicle Routing Problem (SDVRP), along with a few instances, is available as a code template.

Arc Routing

In arc routing problems, the customers to serve are not the nodes of the graph but its edges: a fleet of delivery vehicles with uniform capacity must service edges with known demand. Each edge is served by exactly one vehicle.

The complete Hexaly model for the Capacitated Arc Routing Problem (CARP), along with a few instances, is available as a code template.

Other variants

Many other variants of routing problems (multiple depots, driver breaks, multiple trips, Dial-A-Ride, time-dependent travel times, draft limits, transshipment, …) are available in our routing code templates.