Modeling guide for routing problems

List variables are a very powerful modeling feature for various problems where collections of items have to be ordered in an optimized fashion. For instance, scheduling problems, production planning problems, crew scheduling problems or even assignment problems can be efficiently modeled and solved with list variables.

State-of-the art results can also be obtained for routing problems modeled with list variables. Starting with the traveling salesman problem, this page covers most routing optimization variants and explains how it can easily be modeled with LocalSolver.

Note

The most common routing problems, namely TSP, CVRP and CVRPTW are available in our example tour as complete models with sample data and ready-to-use code in LSP, C++, Java, .NET and Python. This guide refers the reader to these pages and also explores other features like pickup-and-delivery constraints, waiting time minimization, prize-collecting variants, and so on.

The Traveling Salesman Problem (TSP)

../_images/tsp.svg

The simplest routing problem involves only one vehicle, meant to visit a given set of cities with the smallest possible tour.

This problem is easily modeled with a single list variable constrained to contain all the cities. Below is the 3-line LSP model for the traveling salesman problem. Despite its simplicity, this model yields high-quality solutions very quickly. See our TSP benchmark page for a performance comparison with MIP solvers. Complete models for LSP, C++, Java, .NET and Python are available in our TSP example.

// A list variable: cities[i] is the index of the ith city in the tour
cities <- list(nbCities);

// All cities must be visited
constraint count(cities) == nbCities;

// Minimize the total distance
minimize sum(1...nbCities, i => distanceWeight[cities[i-1]][cities[i]])
                + distanceWeight[cities[nbCities-1]][cities[0]];

The Prize-Collecting Traveling Salesman Problem (PCTSP)

../_images/pctsp.svg

Sometimes visiting all the cities is not mandatory. In the prize-collecting traveling salesman problem, a revenue is attached to each city and the objective is to maximize the net revenue defined as the total collected revenue minus the total traveled distance (counting a cost of 1 per unit of distance).

This variant requires only a slight modification of the TSP model, as illustrated below.

// A list variable: cities[i] is the index of the ith city in the tour
cities <- list(nbCities);
revenue <- sum(cities, i => prize[i]);

distance <- count(cities) <= 1 ? 0 : sum(1...count(cities), i => distanceWeight[cities[i-1]][cities[i]])
                + distanceWeight[cities[count(cities)-1]][cities[0]];

maximize revenue - distance;

You can observe that the contraint on the size of the list has been removed. Consequently the distance expression is modified to take into account the now possible case of empty or singleton lists. The revenue collected by the tour is measured with a simple sum. Note that the range on which the sum is applied can be defined as an interval (as in the distance case, to enumerate all positions in the list) or directly as a list (to enumerate values contained in the list).

The Capacitated Vehicle Routing Problem (CVRP)

../_images/vrp.svg

In the capacitated vehicle routing problem (CVRP), a fleet of delivery vehicles with uniform capacity must serve customers with known demand for a single commodity. The vehicles start and end their routes at a common depot and each customer must be served by exactly one vehicle.

The objective is to assign a sequence of customers to each truck of the fleet, minimizing the total distance traveled, such that all customers are served and the total demand served by each truck does not exceed its capacity.

As explained in our example tour, this problem is modeled with one list variable per truck of the fleet, and the coverage of all customers is enforced by a partition constraint:

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

And the capacity contraint is modeled with a sum over all items assigned to a truck:

sequence <- tour[k];
c <- count(sequence);
routeQuantity <- sum(0...c, i => demands[sequence[i]]);
constraint routeQuantity <= truckCapacity;

Many test cases of this capacitated vehicle routing problem are available in the literature. Our CVRP benchmark page shows how quickly LocalSolver finds high-quality solutions to these problems.

The Prize-Collecting Capacitated Vehicle Routing Problem (PCCVRP)

../_images/pccvrp.svg

In the prize-collecting variant of the capacitated vehicle routing problem, instead of serving all customers (which may be impossible due to a limited fleet for instance), the goal is to maximize the total collected revenue.

In this context, the lists will not form a partition of the set of all customers, but it is still forbidden to have more than one truck visiting the same customer. Therefore the partition contraint becomes a disjoint constraint:

// At most one visiting truck per customer
constraint disjoint[k in 0...nbTrucks](customersSequences[k]);

The remainder of the model is unchanged except for the objective function that now measures the total collected revenue, defined on each tour as sum(tour[k], i => prize[i]).

The Capacitated Vehicle Routing Problem with Time-Windows (CVRPTW)

../_images/vrptw.svg

Handling time-windows in routing problems means ensuring that each delivery takes place in a given time-window.

In many cases (for example in the standard CVRPTW), it can be modeled without introducing delivery times as decision variables. Indeed, a key modeling principle is to avoid introducing a decision when an expression can in fact be inferred from other expressions. It is the case for the CVRPTW, where the delivery time for a customer can be defined as the maximum of:

  • the earliest allowed delivery time for this customer

  • the departure time from the previous customer plus the travel time to the current customer

This definition is recursive and can be modeled in LocalSolver by means of a function (i, prev) => .... These functions are used to define the i-th element of an array as a function of the (i-1)-th element and the value i. See our documentation on this topic for details

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

The lateness at each customer is computed as the difference between the ending time of this visit and the latest allowed end for this customer. For the arrival at the depot, we compare the arrival time to the maximum horizon defined for the problem.

Our example tour includes a detailed model for this problem with sample data and ready-to-use code in LSP, C++, Java, .NET and Python.

Multiple time-windows per customer

When dealing with multiple time-windows per customer, we can define a function giving the next available time for visiting the customer.

In the case of two time-windows [A, B] and [C, D], this function taking as parameter the customer index and the arrival time t at customer location can be written as follows:

function nextAvailableTime(customer, t) {
    A <- tw1Start[customer];
    B <- tw1End[customer];
    C <- tw2Start[customer];
    return t <= B ? max(A, t) : max(C, t);
}

The value returned is:

  • A if t is before the first time-window [A, B]

  • t if t is within the time-window [A, B]

  • C if t is between the two time-windows

  • t if t is within the time-window [C, D] or after it (which will lead to lateness in this case)

In the general case with an arbitrary number of time-windows for each customer, the function can be written using the min operator combined with a lambda function:

function nextAvailableTime(customer, t) {
    n <- nbTimeWindows[customer];
    nextTimeWindowStart <- min(0...n, k => twEnd[customer][k] > t ?
            twStart[customer][k] : twStart[customer][n - 1]);
    return max(t, nextTimeWindowStart);
}

This function is used during the computation of end times:

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

Note that you may also precompute an array giving the next available time for each possible time value in order to avoid calling the function directly in the model to improve performance.

The lateness at each customer is computed as the difference between the ending time of the visit and the latest allowed end of last customer time-window.

CVRP with preassignments

../_images/vrp_constrained.svg

Sometimes some visits are preassigned to a given truck. On the contrary, we can also have some visits that cannot be assigned to some trucks (for example because of the size of the truck or because of the clearances of the driver).

Such constraints can easily be defined on list variables thanks to the contains operator:

// with mandatory and forbidden defined as maps of pairs truck/visit
for[req in mandatory] constraint contains(tour[req.truck], req.visit);
for[req in forbidden] constraint !contains(tour[req.truck], req.visit);

The Pickup and Delivery Problem (PDVRP)

../_images/pdp_color.svg

In the pickup-and-delivery problem, the items to be delivered are not loaded into the truck at the depot before starting the tour, but they are collected along the tour. In other words, each demand is made of a pickup point (where the item will be collected) and a delivery point (where the item will be delivered).

This situation induces three constraints:

  • Each pickup/delivery pair must be assigned to the same tour,

  • Within this tour, the pickup cannot be visited after the delivery (an item that is yet to be picked up cannot be delivered),

  • The weight of the truck will vary along the tour (increasing upon pickups and decreasing upon deliveries) and the capacity constraint must be satisfied at any time during the tour.

Considering each list in the LocalSolver model, the first requirement for a pickup/delivery pair (a,b) means that a and b either both belong to the list or none of them belong to the list. Literally it means that contains(sequence[k], a) == contains(sequence[k], b).

As for the ordering of the pair within the list, it is related to the position of a and b within the list, which can be obtained with the indexOf operator:

for[k in 0...nbTrucks] {
    constraint contains(sequence[k], a) == contains(sequence[k], b);
    constraint indexOf(sequence[k], a) <= indexOf(sequence[k], b);
}

Note that indexOf takes the value -1 when searching for an item that is not contained in the list. Consequently writing a strict inequality on the second constraint would be an error because it would not allow both items to be outside the list.

As mentioned above, the last specificity is that the capacity constraint must now be checked after each visit. To do this, we need to define the weight after a visit as the weight before the visit plus the weight of the current item (positive if loaded, negative if delivered):

for[k in 0...nbTrucks] {
    weight[k] <- array(0...c, (i, prev) => prev + demands[sequence[i]]);
    constraint and(0...c, i => weight[k][i] <= truckCapacity);
}

CVRPTW with minimization of waiting time

../_images/vrptw_waiting.svg

In the CVRPTW the usual objective is to minimize the total traveled time. In this context, the satisfaction of time-window constraints is optimal with an earliest scheduling approach (starting the tour at the earliest possible date).

However, this strategy can lead to unnecessary waiting time for drivers. Hence an additional requirement is sometimes to minimize the total waiting time, or (equivalently) to minimize the total duration of tours.

The recommended approach is to introduce a startTime[k] integer decision variable for each tour k. In the CVRPTW model, it just needs to be added to the initial travel time from the depot for each tour:

startTime[k in 0...nbTrucks] <- int(0, maxStart);

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

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

TSP with draft limits (TSPDL)

../_images/draft_limit.svg

The draft of a ship is the distance between the waterline and the bottom of the ship. This draft increases with the load of the ship and each port has a draft limit, beyond which a vessel cannot enter the port. Hence in maritime transportation the ability to visit a site depends on the amount of carried cargo.

Such a limitation can occur in other contexts where the accessibility of a site depends on the weight of the vehicle. As in the pickup and delivery problem, it requires computing the weight of the vehicle along the tour:

weight <- array(0...c, (i, prev) => i == 0
                            ? totalWeight
                            : prev - itemWeight[sequence[i-1]], 0);

// Setting a contraint over all terms of an array is done with the 'and' operator:
constraint and(0...c, i => weight[i] <= weightLimit[sequence[i]]);

// Altenatively (for instance if finding a solution respecting all draft limits
// take more than a few seconds), you can minimize the cumulated overweight:
overweight <- sum(0...c, i => max(0, weight[i] - weightLimit[sequence[i]]));

CVRPTW with a fixed pause

../_images/vrptw_waiting.svg

The CVRPTW with a fixed pause is equivalent to the CVRPTW, but with the addition of a pause of a given length at a given time. The model of the CVRPTW needs to be changed to take this pause into account.

Only a few functions need to be added to take a break during a route. In this problem, during every route, the driver has to take a break of a duration pauseDuration at the exact time pauseTime:

pauseDuration = 60 * 60; // 1 hour break
pauseTime = 12 * 60 * 60; // at midday

It is also considered that if the pauseTime takes place during a service, the pause should be taken before the service. We can then define two functions that compute the different end times for the tasks.

This endTimeSelector function computes the end time of the task depending on its position in the route. It takes into account the earliestStart of the service and the distance to the warehouse to do a pause if necessary:

function endTimeSelector(seq, i, prev) {
    local travelTime <- i == 0 ? distanceWarehouse[seq[0]] : distanceMatrix[seq[i - 1]][seq[i]];
    local lastEndTime <- prev; // prev = 0 if i = 0
    local arrivalTimeWithoutPause <- lastEndTime + travelTime;
    local startTimeWithoutPause <- max(earliestStart[seq[i]], arrivalTimeWithoutPause);
    local shouldBreakBeTaken = lastEndTime < pauseTime
            && startTimeWithoutPause + serviceTime[seq[i]] >= pauseTime;
    local startTime <- max(earliestStart[seq[i]],
           arrivalTimeWithoutPause + shouldBreakBeTaken * pauseDuration);
    return startTime + serviceTime[seq[i]]
}

The returningHomeTime function computes if a pause is needed when the truck is going back to the warehouse by using the endTime of the last task of the route:

function returningHomeTime(last, customer) {
    local returnHomeWithoutPause <- last + distanceWarehouse[customer];
    local isPauseInInterval <- (last <= pauseTime && returnHomeWithoutPause > pauseTime)
            ? pauseDuration
            : 0;
    return returnHomeWithoutPause + isPauseInInterval;
}

We then use these functions in the CVRPTW model and replace the computation of the end of tasks and the homeLateness:

endTime[k] <- array(0...c, (i, prev) => endTimeSelector(sequence, i, prev), 0);

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

The rest of the model is identical to the CVRPTW.

CVRPTW with multiple pauses in time-windows

../_images/vrptw_waiting.svg

The CVRPTW with multiple pauses in time-windows is equivalent to the CVRPTW, but with the addition of multiple pauses that have to take place in a time-window for every truck. The customers also have multiple time-windows to execute the service. The previous model needs to be adapted to take those multiple pauses into account.

Firstly, we need to precise that every pause is of equal duration pauseDuration and that a if a pauseTimes occurs during a service, it should be taken before. The new model needs to decide when to take a break inside the time-window. This requires adding a decision variable for every pause of every truck, in order to decide when the pause starts:

pauseTimes[truck in 0...nbTrucks][p in 0...nbPauses[truck]] <-
        int(startPause[truck][p], endPause[truck][p]);

Secondly, we have to take into account that multiple pauses can occur during a route. In order to do so, we have to write new functions that detect the multiple pauses and add the corresponding duration to the route. The function nextAvailableTime computes the next possible start for a customer by comparing the current time to the start of the time-windows of the customer. It is identical to the function presented in the CVRPTW with multiple time-windows section:

function nextAvailableTime(customer, time) {
    local earliestAvailableTime <- min(0...nbTws[customer], k => twEnd[customer][k] >= time
            ? twStart[customer][k]
            : latestEnd[customer]);
    return max(time, earliestAvailableTime);
}

The function needsPause below simply tests if the pauseTimes is between the beginning and the end of the drive:

function needsPause(pauseStart, start, end) {
    return start <= pauseStart && end > pauseStart;
}

The travelEnd function computes the moment when the travel from one customer to another ends by testing iteratively (more than one pause can occur during the travel) for every pause of the truck if it happens during the drive, and adding the pauseDuration every time:

function travelEnd(vehicle, i, time) {
    local sequence <- customersSequences[vehicle];
    local travelDuration <- i == 0
            ? distanceWarehouse[sequence[0]]
            : distanceMatrix[sequence[i - 1]][sequence[i]];
    local travelEnd <- time + travelDuration;
    local endWithPauses <- travelEnd;
    for [p in 0...nbPauses[vehicle]] {
        endWithPauses <- needsPause(pauseTimes[vehicle][p], time, endWithPauses)
                ? endWithPauses + pauseDuration
                : endWithPauses;
    }
    return endWithPauses;
}

The waitingAndServiceEnd function works the same way to determine the end time of the service by computing if pauses are necessary before and during the service:

function waitingAndServiceEnd(vehicle, customer, time) {
    local nextStartWithoutPause <- nextAvailableTime(customer, time);
    local endWithoutPause <- nextStartWithoutPause + serviceTime[customer];
    local endWithPauses <- endWithoutPause;
    for [p in 0...nbPauses[vehicle]] {
        endWithPauses <- needsPause(pauseTimes[vehicle][p], time, endWithPauses)
                ? nextAvailableTime(customer, pauseTimes[vehicle][p] + pauseDuration)
                + serviceTime[customer]
                : endWithPauses;
    }
    return endWithPauses;
}

The returningHomeTime function is similar to the function for the CVRPTW with a fixed pause returningHomeTime function and is just adapted to take more than one pause if necessary:

function returningHomeTime(vehicle, customer, time) {
    local endWithoutPause <- time + distanceWarehouse[customer];
    local endWithPauses <- endWithoutPause;
    for [p in 0...nbPauses[vehicle]] {
        endWithPauses <- needsPause(pauseTimes[vehicle][p], time, endWithPauses)
                ? endWithPauses + pauseDuration
                : endWithPauses;
    }
    return endWithPauses;
}

Finally, we have to insert the functions inside the model, the same way we did with the CVRPTW with a fixed pause model:

endTime[k] <- array(0...c, (i, prev) =>
        waitingAndServiceEnd(k, sequence[i], travelEnd(k, i, prev)), start[k]);

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

The rest of the model is identical to the CVRPTW.

CVRPTW with regular pauses

../_images/vrptw_waiting.svg

The CVRPTW with regular pauses is similar to the CVRPTW, but a pause has to be taken with at most a certain duration pauseFrequency after the end of the last one. Similarly to the CVRPTW with multiple pauses in time-windows, service can be done in different time-windows.

In order to model this problem, we have to introduce a number of pauses that have to be taken in the period. In order to be sure to have enough pauses to cover the whole period, we can for instance set the number of pauses to:

nbPauses = ceil(maxHorizon / pauseFrequency) + 1;

We can then use decision variables to determine the interval between the end of a pause and the start of the next one. These decision variables can be defined as:

pauseIntervals[k in 0...nbTrucks][p in 0...nbPauses] <- int(1, pauseFrequency);

The pauseTimes defined in the CVRPTW with multiple pauses in time-windows can then be computed by changing them to:

pausesTimes[k in 0...nbTrucks][p in 0...nbPauses] <-
        sum[pause in 0...p + 1](pauseIntervals[k][pause]) + pauseDuration * p;

This allows us to keep the rest of the model almost identical to the CVRPTW with multiple pauses in time-windows. Adding a constraint to force the pauses to cover the whole period is the only required modification. This forbids solutions in which all pauses are at the beginning:

for[k in 0...nbTrucks] constraint pausesTimes[k][nbPauses] >= maxHorizon + 1;

This constraint forces pauses to be taken across the whole period by forcing the last pause to be taken at the end. Other pauses will then be taken throughout the whole period. The rest of the model is strictly identical to the CVRPTW with multiple pauses in time-windows.

Time-dependent CVRPTW

../_images/tdcvrptw.svg

The time-dependent CVRPTW is similar to the CVRPTW, but the travel duration between customers depends on when the travel happens. The horizon is divided in several time ranges. Each time range has its own travel time matrix for travel starting during this range.

To model this problem, we introduce an array timeToMatrixIdx which gives for every time step of the horizon the index of the matrix that should be used .

We also have to define a travel time array with 3 dimensions. The first 2 dimensions represent the origin and destination customers. The third one is the index of the matrix. Hence, travelTime[origin][destination][matrixIdx] is giving the travel time from origin to destination for the matrix at index matrixIdx.

Following the same logic, we define an array that gives the travel time to reach the depot depending on the matrix index. travelTimeWarehouse[customer][matrixIdx] is the travel time between customer and the depot for the matrix at index matrixIdx.

We can then modify the end times defined in the CVRPTW, to select at each step the travel time given by the matrix valid at this time:

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

To compute the home lateness, we slightly modify the definition made in the classical CVRPTW to consider the travel time of the matrix valid after doing the last delivery of the route:

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

Other routing problems

../_images/question_mark.svg

Your vehicle routing problem does not fit in one of the above categories? Do not hesitate to contact us and we will help you with pleasure to find the best model for your problem.