Eﬃcient heuristics for the Rural Postman Problem

A local search framework for the (undirected) Rural Postman Problem (RPP) is presented in this paper. The framework allows local search approaches that have been applied successfully to the well–known Travelling Salesman Problem also to be applied to the RPP. New heuristics for the RPP, based on this framework, are introduced and these are capable of solving signiﬁcantly larger instances of the RPP than have been reported in the literature. Test results are presented for a number of benchmark RPP instances in a bid to compare eﬃciency and solution quality against known methods


Introduction
Consider a weighted graph G = (V, E), with vertex set V = {v 1 , . . ., v p }, edge set E, and edge weights denoted by c(i, j) for all v i v j ∈ E. The well-known Chinese Postman Problem (CPP) is the problem of determining a minimum-weight closed route traversing each edge v i v j ∈ E at least once (Guan, 1962).The CPP is tractable and may be solved in O(|V | 3 ) time (Edmonds and Johnson, 1973).The Rural Postman Problem (RPP) is a generalisation of the CPP in which a subset of the edges E r ⊆ E (called required edges) have to be traversed.It is the problem of determining a minimum-weight closed route traversing each edge in E r at least once.The RPP is NP-Hard (Lenstra and Rinnooy Kan, 1976), except when E r = E, in which case the problem reduces to the CPP.The above CPP and RPP definitions for undirected graphs have been generalised in many ways, and algorithms catering for directed and mixed graphs, for example, have been introduced -see Ball et al. (1995), Dror (2000) and Eiselt et al. (1995aEiselt et al. ( , 1995b) for an overview.These versions of the RPP have many applications -see, for example, Angel et al. (1972), Beltrami and Bodin (1974), Bennet and Gazis (1972), Bodin and Berman (1979), Bodin et al. (1989), Bodin and Kursh (1979), Braca et al. (1993), Desrosiers et al. (1986), Eglese (1994), Eglese and Murdock (1991), Gelders and Cattrysse (1991), Ghiani and Improta (2001), Grötschel et al. (1991), Levy and Bodin (1988), Roy and Rousseau (1989), Stern and Dror (1979) and Wunderlich et al. (1992).
In this paper, local search heuristics are presented for the (undirected) RPP.Heuristic methods stand in contrast to exact methods (i.e., algorithms that yield optimal solutions) in that they are designed with the aim of finding good (though not necessarily optimal) solutions without expending much computational execution time.The class of so-called local search heursitics is a family of methods that operates by iteratively performing transformations (referred to as moves) to existing solutions to an optimisation problem in a way that tends to (but is not guaranteed to) improve the solution as the search progresses.Typically, a local search heuristic operates by considering a number of candidate moves during each iteration, and selects the best one to perform on the solution.During the next iteration, the process is repeated on the transformed solution, and so on.
Perhaps the best known heuristic for the undirected RPP is Frederickson's heuristic (Frederickson, 1979).This heuristic is similar to Christofides' heuristic for the Travelling Salesman Problem1 (TSP), and operates by adding artificial edges (representing shortest paths between the relevant vertices in G) to the subgraph induced by E r in a way that yields a connected, Eulerian graph (i.e., one in which it is possible to find a closed route traversing each edge exactly once).Hertz et al. (1999) Ghiani and Laporte (2000), and Letchford (1996).The algorithm by Ghiani and Laporte has been used to solve instances with 350 vertices to optimality, which seem to constitute the largest instances previously addressed by either an exact or heuristic method in the literature.The heuristic presented in this paper is capable of (approximately) solving considerably larger instances of the RPP, and the quality of its solutions compare very favourably to those documented in the literature for benchmark problems.This paper is structured as follows.A local search framework is described in §2, which allows local search moves that have traditionally been applied successfully to Vertex Routing Problems (VRPs), such as the TSP, also to be applied to Arc Routing Problems (ARPs), such as the RPP and the Capacitated Arc Routing Problem2 .This transition is achieved by the introduction of a complexity reduction method (in §2.3) in which both aspects of the order (in §2.1) and direction (in §2.2) of edge traversals are accommodated.The section is concluded with a worked example, before we report results obtained by the heuristics for a number of benchmark RPP instances in §3.The paper closes with a short conclusion in §4.

Local Search Framework
Denote a solution to the RPP by the sequence S = (v s 1 , v t 1 ), (v s 2 , v t 2 ), . . ., (v sn , v tn ) of required edges in the order in which they are traversed.Traversals taking place between the required traversals are omitted from the sequence and are assumed to take place along routes corresponding to shortest distances between the required edges of the sequence.The total weight of the route is therefore given by where d(k, ℓ) denotes the shortest distance between two vertices v k , v ℓ ∈ V (G), and c(i, j) is the cost weight associated with the edge v i v j , as before.
Typically, one would consider all pairs of these exchanges during a single iteration of the search and then perform one that yields a route of minimum overall cost.In the above example, the traversal order of required edges was altered, but not their traversal directions.However, it may be better to traverse the edge (3,4) in S 1 * , for example, in the direction (4, 3) instead of in the direction (3,4).Consequently, it is necessary to determine the optimal directions of traversals of required edges in S 1 * after performing an exchange.
Applying a move therefore involves altering the order of the required edges in the route, and then determining their directions of traversal.This requirement for determining traversal directions results in an increased time complexity, when compared to applying the same type of move to a VRP.However, by using a complexity reduction method presented later, it is, in fact, possible to determine the cost of a route without redetermining the traversal directions of all of its required edges.This allows for the development of comparatively efficient procedures for many ARPs.
The moves considered in this paper for the RPP are slightly more involved than the above example, and are based on the well-known Two-Opt (Flood, 1956 andCroes, 1958) and Three-Opt (Bock, 1958 andLin, 1965) procedures for the TSP.The application of the Two-Opt method is explained with the aid of Figure 1(a), depicting a solution for a small, hypothetical instance of the RPP.The vertices in the figure represent the required edges of the problem, and the solid arcs represent shortest paths between vertices incident to required edges.The Two-Opt procedure is applied by deleting two of the edges, and reconnecting the route so that it corresponds to the dotted arcs, again using shortest paths between the relevant vertices.Two-Opt is typically implemented as a post-optimisation procedure by performing the best Two-Opt move out of all possible moves, and then repeating this process until no improving solution can be found.
The Three-Opt procedure, depicted in Figure 1(b), is similar to Two-Opt, but three edges are deleted instead and are replaced with shortest paths in a differently connected traversal order.Note that there are eight possible ways of connecting three route segments in this way.Of course the Three-Opt procedure includes Two-Opt as a special case.The vertex 0 in Figure 1 represents a domicile vertex, from which a route is considered to start and end, and may represent a depot in a practical context.However, for the RPP, in which a specific domicile vertex is not specified, the vertex 0 may be taken as any vertex incident to a required edge.
Returning to our previous small example sequence for the RPP, the solution sequence S 2 = (3, 4), (5,4), (5,6), (4, 1), (6,8) is found upon applying the Two-Opt procedure to S 1 , deleting the shortest path between the first and second required edges, and between the fourth and the fifth required edges in S 1 and assuming that the traversal directions of the edges are not also reversed.The Three-Opt procedure for the RPP may be applied in a similar manner.
However, as mentioned earlier, the optimal traversal directions for the edges also have to be determined after performing the standard Two-Opt or Three-Opt operation.

Determining Traversal Directions
The algorithm described in this section may be used to compute optimal traversal directions for the edges in a solution sequence S for the RPP, given a fixed order of the edges in the sequence.The algorithm may be applied after each local search move performed in order to yield the smallest cost weight for the new ordering of S.
Consider a routing sequence S = (v s 1 , v t 1 ), (v s 2 , v t 2 ), . . ., (v sn , v tn ) and construct, from the solution sequence, a directed, layered auxilliary graph L with vertex set The labels b and e represent the vertices at which the route is to begin and end.For a closed route v b = v e , and in the The layered graph L, corresponding to the solution sequence S , respectively] to the edge in L between s i t i and s i+1 t i+1 [t i s i and t i+1 s i+1 , respectively] for every 0 < i < n, where d(i, j) denotes the cost of a shortest path from v i to v j , as before.Similarly assign a weight of d(t i , t i+1 ) [d(s i , s i+1 ), respectively] to the edge in L between s i t i and t i+1 s i+1 [t i s i and s i+1 t i+1 , respectively] for every 0 < i < n.Assign a weight of d(b, s 1 ) [d(b, t 1 ) respectively] to the edge in L between b and s 1 t 1 [t 1 s 1 respectively] and a weight d(t n , e) [d(s n , e) respectively] to the edge between s n t n [t n s n respectively] and e.
Each route from b to e in L represents one way of arranging the traversal directions of edges within S, and a shortest path from b to e in L represents a set of optimal directions by which to traverse the edges of S. For example, if vertex t 2 s 2 is on a calculated shortest path, then the second edge of S should be traversed from v t 2 to v s 2 , and hence coded as (v t 2 , v s 2 ) in S, instead of (v s 2 , v t 2 ).Note that the total weight of a route may be found by adding the sum of the weights of the required edges to the weight of the shortest path.
The computational complexity of finding a shortest path in a directed, acyclic graph, such as L, is O(E(L)+V (L)) (Mehlhorn and Näher, 1999), because no updating of information, such as occurs in Dijkstra's method (Dijkstra, 1959), is necessary during the algorithm execution.However, because no earlier layer of L can be reached from a later layer, and because each layer consists of a predetermined number of vertices, and is connected to the other layers in the particular manner shown, this complexity may be reduced further to O(V (L)).
The algorithm for computing a shortest path between b and e in L is a straightforward one, and is given by the pseudo-code listing in Algorithm 1.Let d L (i, j) be the distance between vertices i and j in the auxilliary graph L, and let w L (i, j) be the weight of the edge between i and j in L. In the algorithm pred(i) is a variable that stores the predecessor vertex of vertex i in a shortest path from b to e in L.
Algorithm 1 (Shortest path through the layered graph L) Input: A layered graph L, as described earlier, with edge weights w L (i, j) for each edge ij ∈ E(L).Outputs: (1) The shortest path from b to e in L, stored in the variables pred(i), for all i ∈ V (L), (2) The shortest distance, d L (b, e), from b to e in L.
2. For i ← 2, . . ., n: Every vertex of L (except b) is considered once by the algorithm, and its complexity is therefore O(V (L)).The operation of assigning the computed distance of a vertex (from some other vertex) is referred to as labelling in the literature on shortest path algorithms.Each step of the algorithm performs a labelling operation on one or two vertices with respect to vertex b.In the next section, it is argued that the above procedure does not necessarily need to be applied each time that a move is evaluated.

Reducing Computational Complexity
The task of computing a shortest path in the layered graph in §2.2 is of linear computational complexity.Hence, if the shortest path is computed each time that a move is evaluated, it adds an order of magnitude to the computational complexity of the search heuristic.However, an approach is presented in this section that allows a move to be evaluated without having to re-compute the shortest path through the layered graph, provided that certain shortest path information is known about the untransformed solution (i.e., the solution before the move was implemented).
For the purposes of describing the method of computational complexity reduction, a move performed in a candidate solution S to obtain a transformed candidate solution S * may be viewed as a sequence of z subsequences, containing pairs, given by where each subsequence represents adjacent elements in S involved in the move.In each pair, o and n represent respectively the old position (in S) and the new position (in S * ) of an edge in S that had its position changed in the sequence during the move.The edges of each subsequence are adjacent both before and after the move (else they are placed in separate subsequences), and the subsequences are listed in increasing order of the positions of their edges after the move, i.e. no subsequence placed later than another will contain a smaller n value than the earlier subsequence, and the pairs within the subsequences themselves are arranged in increasing order of n.To illustrate this notation, consider again our hypothetical solution sequence S 1 in §2.1.Assuming that a move involves moving the second edge to the end of the sequence (which incidentally is one of the moves considered by a Three-Opt procedure), the resulting sequence is (5,6), (5,4), (6,8), (4, 1) , and the move may be expressed as where the second subsequence, consisting of a single pair, represents the edge that is moved from position 2 to position 5, as reflected by its entry (2,5).The first subsequence, consisting of 3 pairs, represents the edges that move to the left to make place for the edge that moved to the end.
that when a subsequence describes edges that are reversed in the transformed sequence, the shortest paths take place along routes that are perhaps not as intuitive as in the case where the edges are not reversed.Consider, for example, the layered graph L The complexity reduction method is now described in general, using the alternative nota-tion for a move introduced earlier in the section.The method is described in some detail in the pseudo-code listing below, and may be read in conjunction with subsequent example, provided after the algorithm listing, to simplify understanding.Define label(L,i,j) as a function that performs the labelling operation of a vertex j (with respect to a vertex i) in order to yield d L (i, j) (see §2.2).Each iteration of a local search procedure proceeds as listed in Algorithm 2.
Algorithm 2 (Local search iteration using the complexity reduction method) Input: A solution sequence S with corresponding layered graph L, as described in §2.2.
Output: A transformed solution sequence, with the minimum total cost out of all the moves considered.
1. Calculate, in L, the shortest distance from all vertices to e (the distances from each of these vertices to all vertices in later layers becomes known).
2. For each possible move (denote the transformed solution S * , and its auxilliary graph L * ): 2.2 For each subsequence k ← 1, . . ., z in TRANSFORM(S → S * ): , more than one pair in subsequence) calculate dL * (b, sn l k tn l k ) and dL * (b, tn l k sn l k ) as follows: (else they are already known) , if the edges of subsequence are not reversed): else (i.e., edges of subsequence are reversed): label(L * , b, sn l k +1tn l k +1), label(L * , b, tn l k +1sn l k +1).
3 Select the move for which dL * (b, e) is a minimum and perform this move on S.
Note that Algorithm 2 does not apply any of the moves during step 2, but computes what the total route cost would be if they were to be applied.The sequence S * is therefore not seen as an existing sequence, but rather as one that may be inferred from TRANSFORM(S → S * ).
The computational complexity of step 1 is O(|V (L)| 2 ).The complexity of step 2 depends on two factors, the first of which is the number of moves evaluated and the second, the time taken to evaluate a move.Using Two-Opt and Three-Opt, for example, there are respectively O(|S| 2 ) and O(|S| 3 ) moves to be evaluated.The method takes O(z) time to evaluate a move, and the complexity of step 3 is O(|S|).Therefore, the worst-case computational complexity of using Two-Opt and Three-Opt in the pseudo-code algorithm shown, is respectively O(|S| 2 ) and O(|S| 3 ).The average-case execution time of the procedures are of the same orders of magnitude.
Finally, although it was stated earlier that TRANSFORM(S → S * ) describes only edges that have their position changed in the transformed solution, in one particular case an edge that does not have its positions changed may be included in the notation to simplify the execution of the complexity reduction method.Consider the Two-Opt move TRANSFORM(S 5 → S 5 * ) = (5, 1), (4, 2) (2, 4), (1,5) performed on a hypothetical solution sequence, S 5 .This move could also be encoded as TRANSFORM(S 5 → S 5 * ) = (5, 1), (4, 2), (3,3), (2,4), (1,5) .The inclusion of the pair (3, 3) does not influence the result of the complexity reduction method, but improves its execution time slightly, because the method evaluates fewer subsequences.Therefore, in a move where a sequence consisting of an odd number of edges is reversed, without otherwise moving their positions in the sequence, it is more efficient to include, in the notation TRANSFORM(S → S * ), the edge that does not have its position changed.Nevertheless, in a computer implementation, this is not a practical concern, because a move would not be encoded as in TRANSFORM(S → S * ) during execution, since this would take O(|S|) time to encode.Rather, the computer program would keep track of the untransformed and transformed positions of the endpoints of each subsequence, and in doing so, would automatically ensure this benefit.The operation of the complexity reduction method is now illustrated by means of example.
The corresponding layered graph, denoted L 0 , is depicted in Figure 6.Its vertex set is given by V (L 0 ) = {1, . . ., 2 × 6 + 2 = 14}.The cost of S 0 may be determined by adding the sum of the cost weights of the required edges (i.e., the value 11) to d L 0 (1, 14) = 15 to yield a total cost of 26.
The evaluation of a single move by the complexity reduction method is described next.The move under consideration is the Two-Opt move TRANSFORM(S 0 → S * 0 ) = (5, 2), (4, 3), (3,4), (2,5) , in which the sequence of edges between the first and the last edges in S 0 are reversed.The resulting layered graph, as it would appear if the move were applied, is shown Figure 7.Note that the edges of this graph are omitted for the sake of clarity.The relevant cost labels required to compute the route weight are shown on the graph, as computed by the complexity reduction method.The computations performed in step 2 of the complexity reduction method proceed as follows: 2 For the only subsequence k ← 1 in TRANSFORM(S 0 ): The total cost of the solution after performing the move would be 24 (i.e., the sum of the cost weights of the required edges, 11, added to d L * 0 (1, 14) = 13).The example illustrates how the total cost of the route may be found without actually performing the move or computing an entire shortest path from vertex 1 to 14.Note also that the optimal traversal directions need not be redetermined when the complexity reduction method is used.4, 3), (3,4), (2,5) .The labels calculated by the complexity reduction procedure are shown as bold numbers beside the relevant vertices.

Computational Results for the RPP
The results obtained from applying the Two-Opt and Three-Opt heuristics (described in §2) to benchmark problems and new test data are presented in this section and contrasted with results previously obtained by others.

Details of Heuristics
A Two-Opt and Three-Opt procedures, described in §2.1, were used to arrive at the results presented in this section.The heuristic of Frederickson (1979), of complexity O(|V | 3 ), was used to generate starting solutions for the Two-Opt and Three-Opt procedures, and the method described in §2.3 was used to reduce the complexity of each iteration of Two-Opt procedure [Three-Opt procedure, respectively] to O(|E r | 2 ) [O(|E r | 3 ), respectively].A one-dimensional array of integer values was used to represent the order in which the edges of the route are traversed, and the Two-Opt and Three-Opt procedures were performed on this array.

Benchmark Graph Instances
The results obtained from applying the Two-Opt and Three-Opt heuristics to a set of benchmark test instances are shown in Table 1.The instances numbered 1 to 24 in Table 1 are the instances used by Christofides et al. (1986), and the other two instances are two of the so-called Albaida instances3 used by Corberán and Sanchis (1991) Note also that instance 18 of the Christofides problems is omitted from this analysis due to inconsistencies in its reported optimal value.Hertz et al. (1999) report that the optimal solution for this graph instance has a cost weight of 147 and that a solution of this cost is found by their heuristics, while Fernández de Córdoba et al. (1998) report that it is 148 and that their procedure also finds a solution of this cost.However, the heuristics described in §2.1 and §2.3 found a solution of cost 146.A manual examination of the route confirms that a feasible solution of cost 146 exists.
Frederickson's heuristic found an optimal solution 9 times out of the 25 instances considered.When the Two-Opt procedure was also applied, this number increased to 15.When the Three-Opt procedure was used in conjunction with Frederickson's heuristic instead, an optimal solution was found 19  ) per iteration) and the fact that actual running times closely match theoretical worst-case times in exchange-based moves (because all exchanges are considered), limits its applicability to small problem instances, in our opinion.

Graph Instances with Euclidean and Random Edge Weights
The graph instances of Table 1 are clearly well within the capabilities of all the heuristics considered, and two larger sets of data were therefore generated, in order to establish of practical feasiblity for the heuristics, in terms of problem instance sizes that could be considered.The results of applying the heuristics to these larger instances are shown in Table 2.
The instances of the first data set in Table 2(a) were generated with Euclidean edge weights 4 , and the instances of the data set in Table 2(b) with random edge weights 5 .The dimensions of the graphs were chosen uniformly in the ranges 500 ≤ |V | ≤ 1 000 and 2 500 ≤ |E| ≤ 3 500, and the number of required edges in each graph (i.e., |E r |) was set to approximately 10% of the size of the graph.The edges designated as required were selected iteratively -a random as yet non-required was chosen and then designated as required only if the subgraph induced by the required edges would be disconnected if it were to be chosen.The process was repeated until the desired number of required edges had been selected.The graphs were generated according to the procedure outlined in Algorithm 3.

Algorithm 3 (Generate Random Graph Instance)
Inputs: Order p and size p − 1 ≤ q ≤ p 2 of problem instance.Output: Random instance of a connected graph G, of order p and size q.
1. Initialise the graph G to contain p vertices (and no edges).
2. Select two different components of G at random, and add an edge between a randomly chosen vertex from each component.Repeat this step until G is connected.
3. While G has fewer than q edges, add an edge between two randomly chosen vertices not yet joined by an edge.
The algorithm used to determine the lower bound values in Table 2 operates in a similar way to the well-known Edmonds & Johnson algorithm for the CPP.A minimum cost maximum cardinality matching is performed on the odd-degree vertices of the subgraph induced by the required edges (Edmonds and Johnson, 1973).
The instances of Table 4 are the largest known RPP instances to which either a heuristic or exact method have been applied successfully.The largest graphs previously reported seem to be those of Ghiani and Laporte (2000), consisting of 350 vertices and an average of 1370 edges (of which approximately 140 are required edges).

Conclusion
heuristics for the Rural Postman Problem (RPP) were introduced in this paper.The heuristics are based on a local search framework that can also be used to optimise routes in other types of arc routing problems (e.g., the Capacitated Arc Routing Problem), and seem to consistute the best heuristic approach currently available for the RPP.

Figure 1 :
Figure 1: Operation of the Two-Opt and Three-Opt move types.Cancelled solid arcs are removed, and the route is reconnected as shown by the dashed lines.

Figure 3 :
Figure 3: The layered graph, L 4 , corresponding to the solution sequence S 4 =

14 Figure 4 :
Figure 4: The layered graph L 4 * of the transformed solution sequence S 4 * .The path indicated by the bold edges, between vertices 5 and 11, represents the same path as the path between vertices 4 and 10 in Figure 3.

Figure 5 :
Figure 5: An example RPP problem instance.Edge cost weights are shown next to each edge, and required edges are denoted in bold face.

times out of 25 .
The heuristic of Fernández de Córdoba et al. (1998) found an optimal solution 11 times out of 25, while the Two-Opt based heuristic of Hertz et al. (1999) found an optimal solution 22 times out of 25.However, the high computational complexity of the Hertz et al. (1999) Two-Opt heuristic (O(|E| 5 introduced a family of local search heuristics for the RPP, while Fernández de Córdoba et al. (1998) employed a heuristic based on Monte Carlo principles.Exact algorithms have been proposed by Christofides et al.

Table 1 :
. Hertz et al. (1999) and Fernández de Córdoba et al. (1998) present results for the same data.Note that the execution times listed in the columns for the Two-Opt and Three-Opt procedures of Table 1 are the total execution time, including the time taken to generate an initial solution.Results for RPP benchmark test instances.Execution times are measured in seconds and were achieved on an Intel Pentium IV processor (2.8 GHz) with 512 Mb memory.