Routing With Sequencing

The motivation for this post comes from a sequence of questions posted on OR Stack Exchange (including this one), having to do with a mixed integer programming model for routing an electronic vehicle (EV) serving various customers. One difference from the basic single vehicle routing models with which I'm familiar is that the EV has to visit a charging station periodically during the route. That is easy to accommodate. Where it gets funky is that the modeler needs to know within the model which customer was last on the route, because the EV is required to go to the nearest charging station after its last stop. I'll take that a step further and require the model to provide (via variables) the position of each node (first, second, third) in the route sequence. This might be useful if, for example, the model had to enforce a rule that customer X must be among the first three customers served. Identifying just the last customer node is easier, as I'll describe at the end.

Attempts by the author of the original question followed the usual pattern for vehicle routing. Assume that there is a single vehicle, each customer must be visited once, and there are no time windows complicating things. You have a digraph containing nodes for each customer and each charging station. You typically start by assigning a binary variable $x_{ij}$ to each arc $(i, j),$ taking value 1 if and only if the vehicle crosses that arc, and proceed from there.

To collect sequencing information, I would normally employ the Miller-Tucker-Zemlin formulation of subtour elimination constraints. The MTZ approach adds a nonnegative auxiliary variable $u_i$ for each node $i$ together with the constraints $$u_j \ge  u_i + x_{ij} - M(1 - x_{ij})$$ for each pair of distinct nodes $i\neq j.$ This says that if we cross arc $(i, j),$ the value of $u_j$ must be at least one higher than the value of $u_i,$ preventing any loops. If $n$ is the number of nodes and we are willing to start numbering with $u_s=0$ ($s$ being the starting node for the tour), we can choose $M=n-1.$ The MTZ constraints are intended to prevent subtours, but as a side effect the $u$ variables number the stops in the order they occur.

This would work for the EV problem if there were a rule that the EV cannot use the same charging station twice during a tour. If the vehicle can stop more than once at the same charging station (which I assume would normally be the case), we cannot use the MTZ constraints because a repeat visit to a charging station would create a subtour.  This also complicates (I think) the use of subtour elimination constraints to prevent disjoint subtours. Fortunately, there are at least two "reasonable" (in my opinion) workarounds, both using the MTZ constraints. Unfortunately both are clunky.

The first workaround is to create multiple clones of each charging station node. So if node $s$ represents a charging station, we introduce additional nodes $s', s'', s''' \dots$ that are all charging stations, all in the same location (meaning time / distance / charge consumption between node $i$ and any of the clones is the same). For any arcs $(i, s)$ and $(s, j)$ we add arcs $(i, s'), (s', j), (i, s''), (s'', j)$ etc. We do not require that every charging node be entered (unlike customer nodes, which must all be visited), but we do limit each charging node to at most one entry. That removes the threat of loops and let us use the MTZ approach. Besides making the digraph larger, this forces the modeler to guess how many clones of each charging node will be needed.

The other approach is change to a multigraph, with fewer nodes but more arcs. We include only customer nodes, plus dummy start and end nodes. Arcs from the start node to each customer (within EV range) are the same as before. The end node is a stand-in for the closest charging station to the last customer visited. The arc from any customer node to the end node has the time / distance / charge consumption required to reach the closest charging station. (The closest charging station to each customer node is computed before building the model.)

For each pair of customer nodes $i \neq j,$ the arc $(i, j)$ (if it exists) represents moving directly from $i$ to $j.$ For select charging nodes $s$ we add another arc from $i$ to $j$, which I will denote $<i, s, j>,$ that represents going from $i$ to $s,$ recharging, and then proceeding from $s$ to $j.$ Each of those arcs produces another MTZ constraint. As usual, every customer node should be entered/exited exactly once.

If the number of charging stations is small, we can create extra arcs  $<i, s, j>$ for every combination of two customers and a charging station, weeding out those that are infeasible (meaning an EV with a full charge could not get from $i$ to $s$ or from $s$ to $j$). To get a smaller model, we can throw out arcs that are dominated by other arcs. For $<i, s, j>$ to dominate $<i, s', j>,$ you would need power consumption from $i$ to $s$ to be no greater than power consumption from $i$ to $s'$ and power consumption from $s$ to $j$ to be no greater than power consumption from $s'$ to $j.$ (If other criteria, such as mileage or transit time, appear in the objective function then they would also factor into the determination of dominance.)

One advantage of this approach is that it would let us dispense with the $u$ variables and the MTZ constraints if the only reason for them was to enable constraints forcing the EV to end at the charging station closest to the last customer, since that is baked into the arcs leading to the dummy end node. If we want to use the cloned charging station approach, we can also use a dummy end node linked to each customer by an arc representing the link to that customer's nearest charging station to enforce the desired ending rule for the tour.

A Bicriterion Movement Model

A question on Operations Research Stack Exchange asks about a bicriterion routing problem. The underlying scenario is as follows. A service area is partitioned into a rectangular grid, with a robot beginning (and eventually ending) in a specified cell. Each move of the robot is to an adjacent cell (up, down, left or right but not diagonal). The robot must eventually visit each cell at least once and return whence it came. One criterion for the solution is the number of movements (equivalently, the amount of time, since we assume constant speed) required for the robot to make its rounds. In addition, each cell is assigned a nonnegative priority (weight), and the second criterion is the sum over all cells of the time of the first visit weighted by the the priority of the cell. In other words, higher priority cells should be visited earlier. Both criteria are to be minimized.

The problem can be modeled as either an integer program or a constraint program. The movement portion is quite straightforward to model. Balancing the two objectives is where things get interesting. One can optimize either criterion after setting a bound on how bad the other can be, or one can use lexicographic ordering of the criteria (optimize the primary objective while finding the best possible value of the secondary objective given that the primary must remain optimal), or one can optimize a weighted combination of the two objectives (and then play with the weights to explore the Pareto frontier). Weighted combinations are a somewhat tricky business when the objective functions being merged are not directly comparable. For instance, merging two cost functions is pretty straightforward (a dollar of cost is a dollar of cost, at least until the accountants get involved). Merging distance traveled and "priority" (or "urgency", or "weighted delay") is much less straightforward. In real life (as opposed to answering questions on ORSE), I would want to sit with the problem owner and explore acceptable tradeoffs. How much longer could a "good" route be if it reduced weighted delays by 1 unit?

I chose to use an integer program (in Java, with CPLEX as the optimization engine), since CPLEX directly supports lexicographic combinations of objectives. You can find the source code in my GitLab repository. A write-up of the model is in a PDF file here, and output for the test case in the ORSE question is in a text file here. The output includes one run where I minimized delay while limiting the distance to a middle value between the minimum possible distance and the distance obtain from lexicographic ordering with delay having the higher priority. It turned out that compromising a little on distance did not help the delay value.

Adding Items to a Sequence

A question posed on OR-Exchange in 2017 asked the following: Given a tour of nodes, how does one best add two new nodes while respecting the ordering of the original tour. Specifically, the author began with a tour 0 - 1 - 2 - 4 - 6 - 0 (where node 0 is a depot) and wanted to add new stops 3 and 5 in such away that, in the revised tour, stop 1 still came before stop 2, stop 2 before stop 4, etc.

This problem can arise not just in vehicle routing but in many sorts of sequencing problems (such as scheduling jobs for production). Of course, preserving the original ordering to the extent possible is not always a concern, but it might be if, for instance, the existing stops are customers who have been promised somewhat general time windows for delivery. In any event, we'll just take the question as a given.

The answer I posted on OR-X made the somewhat charitable (and, in hindsight, unwarranted) assumption that the two new stops would be inserted by breaking two previous arcs, rather than consecutively (for instance, ... - 2 - 3 - 5 - 4 - ...). So I'll post an answer without that assumption here. In fact, I'll post three variants, one specific to the case of adding exactly two stops and the other two more general.

First, let me articulate some common elements. I'll denote the set of original nodes by $N_1$, the set of nodes to be added by $N_2$, and their union by $N=N_1 \cup N_2$. All three approaches will involve setting up integer programming models that will look for the most part like familiar routing models. So we will have binary variables $x_{ij}$ that will take the value 1 if $j$ immediately follows $i$ in the new tour. We will have constraints ensuring that every node is entered and exited exactly once:$$\sum_{j\in N} x_{ij} = 1\quad \forall i\in N\\ \sum_{i \in N} x_{ij} = 1 \quad\forall j\in N.$$The objective function will be some linear combination of the variables (sum of distances covered, sum of travel times, ...), which I will not worry about here, since it is no different from any sequencing model.

The first new wrinkle is that we do not define a variable for every pair of nodes. We create $x_{ij}$ only for the following combinations of subscripts: \begin{align*} i & \in N_{2},j\in N_{2},i\neq j\\ i & \in N_{1},j\in N_{2}\\ i & \in N_{2},j\in N_{1}\\ i & \in N_{1},j\in N_{1},(i,j)\in T \end{align*} where $T$ is the original tour. Thus, for example, we would have $x_{24}$ but not $x_{42}$, nor $x_{26}$. The rationale is straightforward: if we add an arc between two original nodes that were not successors on the original tour, we will force an order reversal. For instance, suppose we replace the arc 2 - 4 with, say, 2 - 6. Node 4 now must appear either before node 2 or after node 6, and either way the order has not been preserved.

Version 1

The first variant makes explicit use of the fact that we have only two new nodes. We add one subtour elimination constraint, to prevent the new nodes from forming a subtour: $x_{35}+x_{53}\le 1.$ Now consider how many different ways we could insert the two new nodes. First, we could break two links in the original tour, inserting 3 in the void where the first link was and 5 in the void where the second link was. Since the original tour had five links there are $\binom{5}{2}=10$ distinct ways to do this. Similarly, we could break two links but insert 5 first and 3 later. There are again ten ways to do it. Finally, we could break one link and insert either 3 - 5 or 5 - 3 into the void. With five choices of the link to break and two possible orders, we get another ten results, for a grand total of 30 possibly new tours.

With that in mind, consider what happens if node 3 is inserted after original node $i$, breaking the link between $i$ and its original successor $j$. (In our model, this corresponds to $x_{i3}=1$.) If this is a single node insertion, then we should have $j$ follow node 3 ($x_{3j}=1$). If it is a double insertion ($i$ - 3 - 5 - $j$), we should have $x_{35}=x_{5j}=1$. We can capture that logic with a pair of constraints for each original arc:
\[ \left.\begin{aligned}x_{i3}-x_{3j} & \le x_{35}\\ x_{i3}-x_{3j} & \le x_{5j} \end{aligned} \right\} \forall(i,j)\in T. \] We could do the same using node 5 in place of node 3, but it is unnecessary. If node 3 is correctly inserted by itself, say between $i$ and $j$, and node 5 is inserted after original node $h$, then the original successor $k$ of $h$ needs a new predecessor. That predecessor cannot be $h$, nor can it be any other original node (given our reduced set of variables), nor can it be node 3 (which now precedes $j$). The only available predecessor is 5, giving us $h$ - 5 - $k$ as expected.

You might wonder how this accommodates a 5 - 3 insertion, say after node $i$. The original successor $j$ of $i$ needs a new predecessor, and 3 is the only eligible choice, so we're good.

I tested this with a small Java program, and it did in fact find all 30 valid revised tours (and no invalid ones).

Version 2

Version 2, which can be applied to scenarios with any number of new nodes, involves building a standard sequencing model with subtour elimination constraints. The only novel element is the reduced set of variables (as described above). A blog is no place to explain sequencing models in their full glory, so I'll just assume that you, the poor suffering reader, already know how.

Version 3

In version 3, we again build a sequencing model with the reduced set of variables, but this time we use the Miller-Tucker-Zemlin method of eliminating subtours rather than adding a gaggle of subtour elimination constraints. The MTZ approach generally results in smaller models (since the number of subtours, and hence the potential number of subtour constraints, grows combinatorially with the number of nodes), but also generally produces weaker relaxations.

The Wikipedia page for the TSP shows the MTZ constraints, although for some reason without labeling them as such. Assume a total of $n$ nodes (with consecutive indices), with node $0$ being the depot. The MTZ approach adds continuous variables $u_i, \,i\in \{1,\dots,n\}$ with bounds $0\le u_i \le n-1$. It also adds the following constraints for all eligible arcs $(i,j)$ with $i\neq 0$:$$u_i - u_j + n x_{ij} \le n-1.$$You can think of the $u_i$ variables as counters. The MTZ constraints say that if we go from any node $i$ (other than the depot) to any node $j$ (including the depot), the count at node $j$ has to be at least one larger than the count at node $i$. These constraints preclude any subtours, since a subtour (one starting and ending any place other than the depot) would result in the count at the first node of the subtour being larger than itself.

As I mentioned, the MTZ formulation has a somewhat weaker LP relaxation than a formulation with explicit subtour elimination constraints, so it is not favored by everyone. In our particular circumstance, however, it has an additional virtue: it gives us a relatively painless way to enforce the order preservation requirement. All we need do is insert constraints of the form$$u_j \ge u_i + 1\quad\forall (i,j)\in T.$$This forces the counts at the original nodes to increase monotonically with the original tour order, without directly impacting the counts at the new nodes.