Solving the hierarchical windy postman problem with variable service costs using a math-heuristic algorithm

The Hierarchical Windy Postman Problem (HWPP) is an arc routing problem in which an order relation is imposed on the arcs/edges of the graph, and one has to pass through each edge at least once while adhering to the hierarchical priority relations. The tour starts from and ends at a specific node and the aim is to minimize the length of the tour. We consider a variant of the HWPP in which (i) the precedence order of the edge hierarchies is linear and edges within each hierarchy are connected and (ii) the cost of serving each edge decreases with the number of times it is traversed, and we refer to it as HWPP with variable service costs. An integer non-heuristic linear mathematical formulation is proposed, and a solution approach is designed. Our solution heuristic adapts the layer algorithm of Dror et al. (Networks 17:283–294, 1987) but employs an integer mathematical formulation as a sub-procedure instead of the blossom algorithm to find the least cost path between the nodes of the graph. This choice is based on the fact that the blossom algorithm requires a symmetric cost structure while we deal here with the general case of asymmetric cost structure, which makes our problem a windy variant of the postman problem. It should be noted that our problem is not asymmetric in the sense that there are no opposite arcs with different costs but there are edges which have different costs depending on the traversal direction. In order to compare the performance of our heuristic algorithm with respect to the performance of the mathematical model that is solved by the commercial solver Gurobi, 84 test instances are generated having varying sizes and densities and with different number of hierarchies. These test instances are solved by both methods and the generated results show that the proposed heuristic method is much faster and generates better quality solutions.


Introduction and literature review
Arc Routing Problems (ARPs) are a broad class of problems in the graph theory that arises in delivering many kinds of services. The recent surveys of Lahyani et al. (2015), Mourão and Pinto (2017) and Corberán et al. (2021) enumerate several applications of ARP including street cleaning, salt spreading, snow cleaning, road maintenance and drone inspection. If capacity limitations are imposed on the arcs, the problem is called as capacitated ARP (see Golden and Wong, (1981) for the first introduction of capacitated ARP or see Mandal et al. (2015) for a more generalized version of capacitated ARP). If services of the arcs have to follow certain periodic patterns, the problem is called as periodic ARP (see for instance Monroy et al. (2013) and Triki (2017)). Mathematically, ARP can be defined on a connected graph G ¼ ðN; EÞ having n nodes and m arcs/edges and where N is the node set and E represents the arc/edge set. An edge set can be said to exist if traversal is possible in both directions, otherwise, if transition is only possible in the intended direction, then there is an arc set at hand. The graph can be either directed or undirected and to every arc/edge is associated a cost that can represent, for example, its length or traversal time. The focus of ARPs is on servicing the arcs/edges of the graph, rather than serving its nodes. A particular case of ARP, known as the Chinese Postman Problem (CPP), consists in identifying a least cost closed tour that starts from a depot, traverses each arc/edge of the graph at least once before coming back to the same depot. The CPP is usually defined over an undirected graph and all its edges need to be serviced. Some of the abovementioned applications may be characterized by the fact that the edges do not have the same importance of being serviced but are rather categorized according to their class of priority (called hierarchies). We deal here with the Hierarchical Windy Postman Problem (HWPP) that has the same features of the HCPP (Hierarchical Chinese Postman Problem) with the only exception that the service costs do not have to be symmetric (see Keskin and Triki 2022).
The HCPP is defined in Afanasev et al. (2021) as the problem of finding the shortest traversal of all edges of a graph, respecting the precedence constraints given by a partial order on classes of edges. The authors also identify three special variants of the problem: (i) the case in which the precedence order of the edge hierarchies is linear, denoted as HCPP(l) (ii) the case involving connected edges within each hierarchy, denoted as HCPP(c) and (iii) when both (i) and (ii) restrictions hold, and was denoted as HCPP (c,l). In the sequel we will refer to third variant of the problem, i.e., HCPP(c,l) where a linear precedence order is imposed and each edge hierarchy induces a connected subgraph, unless differently specified. Likewise, we will often refer to the HCPP terminology in the problem definition and the literature review sections since most of the features that ground on HCPP apply to the HWPP as well.
It is to be noted that especially the second restriction related to the connectivity of the edges within each hierarchy does not to seem too restrictive from the practical viewpoint. Indeed, in the urban snowplow problem, for example, interconnected main roads have primary priority to be plowed which will form the first hierarchy. Similarly, second hierarchy consists of second priority roads that are usually connected and so on. Thus, our assumption that each hierarchy is connected is fully compatible with reallife applications.
Formally, if we denote by jHj the number of pre-defined hierarchies and by E 1 ; E 2 ; . . .; E jHj the subsets of edges belonging to each hierarchy, then all edges belonging to E g must be serviced before any one belonging to E h (with E g T E h ¼ £) if g has a higher hierarchy than h, i.e., if h; g 2 f1; . . .; jHjg and g\h. However, there are contradicting views among the scholars on the feasibility of traversing any edge (with or without servicing it) having a lower hierarchy before completing the service of any other edge belonging to a higher hierarchy category. This fact will define, indeed, three different variants of the HCPP: • The first allows any edge from class E h to be traversed (without service) even before completing the service of the edges belonging to a higher priority class E g , i.e., when g\h. This variant of the problem, that can be applied for example to the case of street cleaning and inspection, meter reading, waste collection, and restoration of the roads after disasters, has been proposed by Alfa and Liu (1988), Eiselt et al. (1995), Colombi et al. (2017a), Colombi et al. (2017b), Shao et al. (2020), Oruc and Kara (2018) and Akbari et al. (2021) and is, sometimes, called HCPP with weak precedence relation. • The second variant of the HCPP adopts a more restrictive interpretation of the precedence relation (strong relation) and does not allow any edge from class E h to (not even) be traversed before completing the service of all edges of class E g having a higher priority. This interpretation of the hierarchy has been adopted by Dror et al. (1987), Ghiani and Improta (2000) and more recently by Ç odur and Yılmaz (2020) and can be suitable, for example, for the snow removal application since it would be difficult traversing a road that still needs to be cleaned. • Yet, a different interpretation of the precedence relation has been introduced by Quirion-Blais et al. (2017) and that arose in a case-study related to the multi-vehicle road snow plowing and de-icer spreading. It consists in allowing the vehicles to start plowing any category of streets but those having higher priority should be completed before all the others.
In this paper, we will adopt the second and more restrictive type of precedence relation, but our mathematical model as well as our heuristic method can be easily adapted to suit the weak hierarchy requirement. We will consider that the cost of serving each edge will decrease with the number of times it is traversed (see Corberán et al. (2014) and Keskin and Yılmaz (2019)). This variant of the problem, that we call the Variable HWPP (V-HWPP), belongs to a wider family of variable service cost routing problems. In the sequel, we will first review the most relevant works that dealt with the different variants of the HCPP and the approaches adopted for their solution, and then we will analyze the different ways of assuming variable service cost and identify the specific version we selected for our study.

Related works on the HCPP
The HCPP has been defined over three decades ago in the seminal paper of Dror et al. (1987). The authors dealt with a CPP in which the edges of the underlying network are roads characterized by a hierarchical precedence relation.
Even though their HCPP defined on an undirected graph is NP-hard, they succeeded to solve for the special class HCPP(c,l) in O(kn 5 ), where n is the cardinality of the nodes set and k is the number of hierarchies in the network. Later on, several studies have focused on reducing the complexity of Dror et al. (1987)'s algorithm by trying to improve its features. Ghiani and Improta (2000) developed an exact approach for the HCPP(c,l) having complexity O(k 3 n 3 ) by constructing an alternative auxiliary graph. In this same direction, Korteweg and Volgenant (2006) reduced the complexity of Dror et al. (1987)'s algorithm to O(kn 4 ) but rather than focusing on building a different auxiliary graph, they employed post-optimality techniques to improve the matching in the resulting auxiliary graph. Given the complexities of both the above exact approaches, the selection of the most appropriate algorithm depends on the characteristics of the instance to be solved, i.e., the value assumed by the number of nodes n compared to that of number of hierarchies k (Korteweg 2002).
Similar to Ghiani and Improta (2000), Cabral et al. (2004) developed an exact solution approach but for the HCPP(l) that was based on the idea of transforming the HCPP into a Rural Postman Problem (RPP) and solving it by a branch-and-cut technique (see also Ghiani et al. 2005Ghiani et al. , 2008. The attention of the researchers in solving the HCPP has been directed even toward developing heuristic approaches. The first heuristic method, proposed by Alfa and Liu (1988), was defined for a directed network with weak precedence relation and without any connectivity restrictions on the arcs belonging to the same hierarchy. Their method consists in connecting and balancing the nodes that are endpoints of arcs belonging to the same hierarchy and then producing a single CPP tour based on the hierarchywise tours. Later on, Alfa and Liu's method was outperformed by another heuristic suggested by Damodaran et al. (2008) which was based on a more efficient way of combining the unconnected subnetworks. Damodaran et al. (2008) focused on proposing good quality lower bounds for the HCPP with the aim of defining an efficient heuristic method. Colombi et al. (2017a) employed a math-heuristic to solve the HCPP starting from the exact solution of a Mixed RPP defined on an auxiliary graph for each hierarchy, whereas Colombi et al. (2017b) performed a polyhedral analysis of a hierarchical mixed rural postman problem and developed a branch-and-cut algorithm for its exact solution. For the sake of completeness, it is worth noting that meta-heuristic approaches, namely a genetic algorithm and a hybrid simulated annealing, have also been employed by Ç odur and Yılmaz (2020).
The last family of approaches defined to solve the HCPP is based on solving the minimum spanning tree problem. Sayata and Desai (2015) implemented an efficient heuristic based on the use of Kruskal's algorithm to reduce the number of edges in the underlying auxiliary graph. More recently, Colombi et al. (2017a) tackled the HCPP with mixed graph and employed a minimum spanning tree problem to identify an initial feasible solution and then embedded it within a tabu search procedure in order to improve and diversify the generated solutions.

Variable service cost in HWPPs
As mentioned above, one of the novel aspects of our work is related to considering a particular variant of the HWPP in which a variable service cost scheme is implemented, i.e., the V-HWPP. Most of the scientific literature in the context of routing optimization has focused on variable service cost expressed in terms of time varying models. Interested readers may refer to the excellent review of Gendreau et al. (2015), in which the authors highlighted the scarcity of contributions related to the time varying ARPs, compared to their node routing counterpart. Few works appeared after Gendreau et al.'s review to deal with the time-dependent ARP in general (such as Sun et al. (2015) and Vidal et al. (2021)) and one is specifically dealing with the time-dependent HCPP (Ç odur and Yılmaz 2020). However, our contribution here is not related to the variable service cost from the traversal time viewpoint, but rather to the service cost that decreases with the number of server passages. To the best of our knowledge there are only two papers that have adopted this particular variable cost scheme in the context of the windy CPP. The first is due to Dussault et al. (2013) who solved a snow plowing problem in which the cost of traversing a street varies based on its order within the plowing tour. Specifically, the authors considered that the cost of serving the roads decreases after the first pass by the plowing vehicle. The second paper, by Keskin and Yılmaz (2019), generalizes variable service cost scheme by assuming that each server pass will contribute to reducing further the roads' service cost. This assumption seems to be very realistic and intuitive in a wide range of applications in which the multipassing through roads and/or the learning process by the server contribute to decreasing the deadhead cost. For instance, the postman is expected to learn the roads better as the number of passes increases in the post-delivery application. Similarly, if each pass on the edges represents the shoveling of the accumulated snow on a specific edge, then the cost of the shoveling should not increase by the number of the passes on that edge as the snow becomes less and less at each pass of the snowplow vehicle. Our study will, thus, adopt the variable cost scheme of Keskin and Yılmaz (2019) and will implement it for the first time in the framework of the HWPP. On the other hand, it is a common practice to differentiate the service cost of the arcs from deadheading cost of the arcs as it is done for instance in Perrier et al. (2008). Similarly, Aráoz et al. (2013) associate prizes to arcs rather than costs and also differentiate between service and deadheading prizes, namely they take deadheading prize as zero.

Problem definition and our contributions
This study will specifically focus on solving the V-HWPP, a particular variant of the HWPP in which is defined over an asymmetric network whose service cost decreases with the number of edge passes. Moreover, we restrict our attention here to the variant of the problem characterized by a linear precedence order of hierarchies and by a connected subgraph within each hierarchy.
From the above analyses it results evident that the above V-HWPP has never been addressed in the scientific literature. The objective of this study is to fill in this gap and to provide the practitioners and academic researchers with a new methodology that can contribute to solving the several real-life applications that arise in this specific ARP context.
The paper is organized as follows. Section 2 will formally introduce the problem and will propose an original formulation for the V-HWPP. Then an ad hoc mathheuristic approach based on the concept of layer algorithm and computational experiments are discussed in Sects. 3 and 4, respectively. Finally, some concluding remarks and future research suggestions are provided in Sect. 5.

Optimization model
In this section, we first describe the sets, parameters, and variables that are used in the formulation. Later, the mathematical formulation of the V-HWPP is provided.

Sets, parameters and variables
We assume that the set of nodes and the set of edges connecting the nodes are, respectively, denoted by N and E. Indices i and j are used to represent nodes of N and an edge connecting nodes i and j is denoted by edge ði; jÞ. H represents the set of hierarchies and index h is used to specify a hierarchy. The set of edges that belongs to hierarchy h is denoted as E h . T stands for the set of steps while indices t and m are used to refer to specific steps while jTj represents the maximum number of steps. By the way, a step is used as a measurement tool to track the transition of the vehicle from one node to another. In other words, at each step, the vehicle passes from one node to another. For instance, if a route of a vehicle consists of nodes 1, 5, 3, 2, 5, 1, then we say that the vehicle is at the starting node 1 at the beginning of step 1 and it goes from node 1 to node 5 at step 1, it is at node 5 at the beginning of step 2 and it goes from node 5 to node 3 at step 2, then it is at node 3 at the beginning of step 3 and it goes from node 3 to node 2 at step 3, etc. Note that, the duration of the steps may change from step to another and from a vehicle to another as the distances between the node-pairs do not have to be identical. Finally, K is the set of number of passes and the index k is used to point out the members of K.
The service cost parameters used in the formulation depend on the above-described sets and indices. c ijk denotes the k th service cost of edge ði; jÞ if the service is through the direction from i to j. Note that we assume an undirected graph implying that edges can be traversed through either direction. However, if an edge is served (independently from the service direction), the cost of the service reduces for both directions which also makes sense for real-life applications. For instance, if a road is served for a snowplow operation then the cost of the next service reduces for both directions. In addition, we assume that c ijk 1 ! c ijk 2 if k 1 \k 2 . That is, service costs do not become costlier as the number of the passes increases. Note that for each edge ði; jÞ we also define c jik in the same manner to define the service cost for the direction from j to i. Finally, parameter jE h j represents the number of edges belonging to set E h .
There are five sets of decision variables in our model. h ht , x ijt , y ijt , and a ijtk are the binary variables of the model. h ht shows whether all edges of hierarchy h are served before step t, or not. x ijt (similarly x jit ) indicates whether or not edge i; j ð Þ is served at step t through direction from i to j (through direction from j to i), and y ijt takes value 1 if edge ði; jÞ is served before or exactly at step t independently from the passage direction, and it is 0 otherwise. Finally, a ijtk is 1 if edge ði; jÞ is served exactly k times before or at step t independently from the direction, and it is 0 otherwise. Our model results to be an integer nonlinear program since the objective function involves a multiplication of the binary variables a ijtk and x ijt (a ijtk and x jit ). For the sake of convenience, a summary of all the sets, parameters and variables used along our formulation are provided in Table 1.

Mathematical Model
Our integer nonlinear program abbreviated as V-HWPPF (Variable HWPP Formulation) is given as: The objective function (1) minimizes the total service cost. Constraint (2) and constraint (3) ensure that vehicle leaves the first node (depot) at step 1 and goes back to the same node at the last step, jTj. Note that jTj is an upper limit value for the number of steps that can be considered in any day (we selected its value as the number of edges times twice the number of hierarchies). In order to avoid forcing the vehicle to perform exactly jTj steps on each day, we added a dummy arc that connects the starting node to itself and that has a traversal cost equal to 0. Thus, if on any day the vehicle needs to return to the starting point with fewer steps than jTj, it can return to the starting node at the end of step jTj by traversing the dummy arc at 0 cost during the remaining steps. This ensures the satisfaction of constraint (3) for that day. Constraint (4) guarantees that the vehicle travels through a single edge at each step. By constraint (5) we make sure that each edge is traversed at least once independently from the direction of the passing. Constraint (6) ensures that the vehicle leaves any node i at step t þ 1 if it entered to it at the previous step t, and moreover it cannot leave node i at step t þ 1 if it did not enter it at the previous step t. Constraint (7) determines the value of each variable a ijtk that takes value 1 depending on the number of travels through each edge before or at step t. The constraint ensures that if a ijtk is 1, then the number of travels through the edge (independent of the direction of the passing) is equal to k. Note that only one of the a ijtk variables existing on the left-hand side of the constraint will be 1 in the optimal solution, although it is possible to satisfy the constraint by setting more than one a ijtk variables to 1. This is because the objective is a minimization Indicates whether or not all edges belonging to E h are served until step t x ijt (x jit ) Indicates whether or not edge ði; jÞ is served at step t through the direction from i to j (through the direction from j to i) y ijt Indicates whether or not edge ði; jÞ is served before or at step t in total, independently from the direction a ijtk Indicates whether or not edge ði; jÞ is served k times before or at step t independently from the direction Solving the hierarchical windy postman problem with variable service costs using a math-heuristic… 8793 function and the variable service costs are chosen to decrease by each pass. For instance, suppose that the number of total passes of an edge i; j ð Þ 2 E before or step t is 5. Hence, the right-hand side of constraint (7) written for edge i; j ð Þ and for period t will be 5. It is possible to assign, for instance, both a ijt2 and a ijt3 to 1 (and all other a ijtk variables to 0) which will make the left-hand side of the constraint as 2 Â a ijt2 þ 3 Â a ijt3 which is equal to 5. Another alternative is, of course, to set only a ijt5 variable to 1, making the left-hand side of the constraint as 5 Â a ijt5 which is also equal to 5. Since the cost of the latter alternative in the objective function, which is c ij5 , is lower than the cost of the first alternative, which is c ij2 þ c ij3 , the solver will naturally choose the latter alternative. Note that if an edge, say i; j ð Þ; is traversed through i to j, then the next passage cost should be decreased not only for the passage direction from i to j but also for the passage direction from j to i. Hence, if an edge i; j ð Þ 2 E is traversed for instance 4 times in total independent of the passage direction, its fifth passage cost will be c ij5 if the passage direction is from i to j and c ji5 if the passage direction is from j to i. Constraints (8) and (9) define the relationship between x ijt , x jit variables and y ijt variable. Namely, if the number of travels through an edge is zero until a particular step, i.e., [ 0 must hold. Constraint (10) defines the variable h ht , i.e., it guarantees that h ht variable is equated to 0 if if any edge belonging to hierarchy h is not traversed before or during step t. Similarly, constraint (11) avoids traversing edges belonging to hierarchy h þ 1 if all the edges of hierarchy h are not traversed yet. Finally, constraint (12) puts the usual binary restrictions.
The above nonlinear model can be easily linearized by introducing a new continuous variable l ijkt (l jikt ) to replace each multiplication a ijtk x ijt (a ijtk x jit ). We also introduce seven new sets of constraints, which are linear and force l ijtk (l jitk ) to be equal to a ijtk x ijt (a ijtk x jit ). These constraints are given in the following.
Note that if there is only one hierarchy and if the cost parameters are taken as constant then the V-HWPP reduces to the WPP which is known to be an NP-Complete problem. Hence, V-HWPP is also a difficult problem and the solution time for formulation (1)-(19) will increase exponentially with the size of the instances which is characterized by the number of nodes, edges and hierarchies. For this reason, we will propose in the next section a heuristic method that is able to tackle even large-scale instances.

Windy layer algorithm heuristic
In our heuristic approach, we adapt the layer algorithm of Dror et al. (1987), that solves the HCPP with linear precedence relationships and edges connected within each hierarchy to our variant of the problem. Note that the exact solution methods of Ghiani and Improta (2000) and Korteweg and Volgenant (2006) are more efficient in terms of computation time, but their methods also depend on the layer algorithm as well. The layer algorithm employs the famous blossom method of Edmonds (1965) as a subprocedure to find the path that goes through all the edges of the current hierarchy without passing through any edge from the lower hierarchies. In the sequel, we will call, for convenience, this path as the least cost path. However, the blossom algorithm requires the edge traversal costs to be symmetric and as our problem is a windy postman problem the edge traversal costs possess an asymmetric nature. Hence, we construct here mathematical models and run them as sub-procedures in order to find the least cost paths, passing through each edge belonging to the current hierarchy, between any two nodes within the layer algorithm. We will refer to this heuristic as the Windy Layer Algorithm (WLA). The details of such least cost path mathematical models are presented in the following after explaining the WLA in the context of the V-HWPP.
Let G h ¼ ðN h ; E h Þ be the subgraph induced by the set of edges of hierarchy h. Here, N h represents the incident nodes of E h . Define F q as the subgraph induced by the first q sets of edges, i.e., F q ¼ S q h¼1 G h . We assume that subgraphs F q are connected for all q ¼ 1; . . .; jHj. Moreover, we define node i 2 N h as an entry node of G h if it is also incident to an edge of F hÀ1 . The set of entry nodes of G h is Finally, we define Y 1 ¼ fi 1 g and Y jHjþ1 ¼ fi 1 g where i 1 is the depot which is the initial and the final node. Now, we construct a new graph h¼1 Y h . Note that if a node exists in more than one of the sets Y h ; h ¼ 1; . . .; H j j þ 1, it will be treated as a different node each time. Hence, the nodes will exist exactly the number of times they exist in the sets Y h ; h ¼ 1; . . .; H j j þ 1. On the other hand, the edge set is defined as Namely, we define an edge for each pair that belongs to successive Y h sets. The graph G 0 is depicted in Fig. 1.
Note that if edge i; j ð Þ 2 E 0 then i 2 Y h and j 2 N hþ1 for some h ¼ 1; . . .; H j j and going from i to j in G 0 means that each edge belonging to E h is traversed starting from node i and ending at node j by following the least cost path (edge lengths are taken as service costs) without violating the hierarchy restrictions. Therefore, the cost of each edge in E 0 , i.e., cost of edge i; j ð Þ, is determined by finding the path having the least cost between node i and node j while traveling through every edge of E h without passing through any edge from lower hierarchies. As mentioned earlier, instead of employing the blossom algorithm, we make use of mathematical programs in order to find the least cost paths between the nodes. Let Pði; j; hÞ denotes the problem of finding the least cost path between node i and node j that travels through every edge of E h without passing through any edge from lower hierarchies. Now suppose that u and v represent the nodes for which u; v ð Þ 2 E and b i uv (b i vu ) indicates the number of times edge u; v ð Þ is traversed through direction from u to v (from v to u) along the least cost path from the starting point until node i without violating the hierarchy restrictions. Note that Pði; j; hÞ is an open WPP and any edge may be traversed more than twice in the optimal solution. Nevertheless, we restrict the number of traversals to two in order to obtain a very fast solution of Pði; j; hÞ since this problem will be solved very often throughout the algorithm (i.e., at every step corresponding to different values of i, j and h). We are aware that it is possible, for some instances, that a vehicle may need to pass some edges more than twice in the optimal solution. However, we observed through preliminary empirical experiments that the rate of such occurrence is very low and restricting the number of passes to two will not harm the optimality for most of the instances while contributing in significatively reducing the overall solution time. Moreover, we drop index t from x ijt variables in the problem. Consequently, we define x uv (x vu ) as a continuous variable that indicates the number of times edge u; v ð Þ is traversed through the direction from u to v (from v to u). We also define two new sets of variables y uv and z uv (y vu and z vu ) that, respectively, indicate if u; v ð Þ is traversed once or twice through the direction from u to v (from v to u) in each solution of Pði; j; hÞ. The formulation of problem Pði; j; hÞ with these new set of variables, for each i, j and h, is given below.
s.t. Solving the hierarchical windy postman problem with variable service costs using a math-heuristic… 8795 x vu ¼ y vu þ 2z vu u; v ð Þ 2 E ð26Þ x vu 2 u; v ð Þ 2 E ð30Þ We minimize the total traversal cost in (20). Note that an edge belonging to higher hierarchies (i.e., an edge from . . .; h À 1), may have been passed several times before we start to traverse the edges of E h : As the traversal costs of the edges vary with the number of passes through edges, we keep track of the number of passes for each edge and for each direction. If edge ðu; vÞ is traversed through u to v b i uv times before, its next traversal will be the (b i uv ?1) st time. Hence, y uv variable (which indicates whether or not edge ðu; vÞ is traversed from u to v once in Pði; j; hÞ, i.e., during traversing edges of E h without passing through any edge from lower hierarchies starting from node i and ending with node j) is multiplied by the corresponding cost c uvðb i uv þ1Þ . Similarly, z uv variable (which indicates whether or not edge ðu; vÞ is traversed from u to v twice in Pði; j; hÞ) is multiplied by the corresponding total cost which is ðc uvðb i uv þ1Þ þ c uvðb i uv þ2Þ Þ. In a similar manner, the objective function coefficients related to y vu and z vu variables are, respectively, determined as c vuðb i vu þ1Þ and ðc vuðb i vu þ1Þ þ c vuðb i vu þ2Þ Þ. However, we still need to adjust the objective function value after obtaining the solution. The point is that after passing through an edge ðu; vÞ from u to v, the cost of passage from v to u should also be decreased and vice versa, and this reality is not captured in the objective function (20). Hence, after solving Pði; j; hÞ problem, we determine the exact route of the vehicle making use of the values of the x, y and z variables and adjust the objective function value accordingly. Therefore, after solving Pði; j; hÞ for each edge i; j ð Þ in E 0 for each hierarchy h, we adjust the objective function value depending on the route of the vehicle obtained from the x uv and x vu values. One may object here that Pði; j; hÞ should have been constructed in such a way that the objective function is calculated correctly. However, in order to be able to drop the t index from x variables and in order to construct Pði; j; hÞ without using the a variables (that has four indices) we had to construct Pði; j; hÞ formulation as given by (21)-(33) and adjust the objective function later. This is a reason for which our method is not an exact method but a heuristic. By constraint (21) we make sure that each edge from E h is traversed at least once independently from the direction of the passing. Constraint (22) maintains the flow balance for each node different from i and j. In other words, constraint (22) guarantees that the total number of enters to a node different than i and j is equal to number of exits from the same node. On the other hand, as the journey starts from node i and finishes at node j for Pði; j; hÞ, the total number of exits minus total number of enters should be equal to 1 and -1 for nodes i and j, respectively. These requirements are satisfied by the help of constraints (23) and (24), respectively. Constraint (25) and (26) define y uv and z uv (y vu , z vu ) variables and sets their relationship with x uv (x vu ) variables. Constraints (27) and (28) avoid passing through an edge belonging to lowerlevel hierarchies and constraints (29) and (30) limit the number of passages of edges by 2 for both directions. Finally, constraint (31) defines x uv and x vu as continuous variables while constraint (32) defines y uv , y vu , z uv , and z vu as binary variables. It should be noted that each x uv and x vu can take only the integer values 0, 1 or 2, but we can define them as a continuous variable (to improve the efficiency of the model solution) knowing that they cannot take fractional values due to constraints (25) and (26). Note that although problem Pði; j; hÞ is an open WPP, it does not require subtour elimination constraints. This is because the edges of each hierarchy are connected within themselves. On the contrary, suppose there is a path from i to j including some (or all) edges of hierarchy h and probably some edges from higher hierarchies and a subtour that is disconnected from the path. If the disconnected subtour does not include any edge from hierarchy h, then it will not be included in the optimal solution since it is not necessary to pass through edges included in the subtour. Hence, any subtour that does not include any edge from hierarchy h will be eliminated by the solver. On the other hand, if there is at least one edge from hierarchy h in the subtour, then there must be at least one edge from hierarchy h connecting the subtour with the path that makes the subtour and the path connected which is a contradiction. This implies that although formulation (20)-(32) does not include subtour elimination constraints, disconnected subtours are not formed in the optimal solution due to the connected nature of the edges within the hierarchies.
By running Pði; j; hÞ for each i; j ð Þ 2 E 0 such that i 2 Y h and j 2 N hþ1 for each h ¼ 1; . . .; H j j, the length of every edge of graph G 0 can be calculated, i.e., the length of each i; j ð Þ 2 E 0 such that i 2 Y h and j 2 N hþ1 for some h ¼ 1; . . .; H j j is the optimal value of the related Pði; j; hÞ. Then, a shortest path from i 1 2 Y 1 to i 1 2 Y H j jþ1 can be easily found by Dijkstra's label algorithm (Dijkstra 1959). We do not think it is necessary to give details of this very wellknown algorithm here. However, along the algorithm, it is necessary to calculate the number of times each edge is passed through each direction while proceeding on the shortest path from the beginning node to each point. In order to achieve this, if shortest path (from i 1 2 Y 1 to i 1 2 Y H j jþ1 in graph G 0 ) pass through node i just before node j, where x uv and x vu values come from the solution of P i; j; h ð Þ. By this updating mechanism, it is ensured that each b i uv is always equal to b i vu for every edge u; v ð Þ and for each i 2 N.
The layer algorithm in which problem P i; j; h ð Þis solved by the above-mentioned mathematical model could have been an exact solution algorithm for the V-HWPP. However, the solution of the P i; j; h ð Þ problems take more and more computational time as the size of the original network increases. Hence, during the implementation phase, we avoided traversing the edges more than twice (as mentioned before) and we forced the commercial solver Gurobi to finalize the run as soon as it finds a feasible solution having an objective function value that is at most 5% away from the optimal value and returns the feasible solution it finds. Consequently, the specific variant of the layer algorithm WLA that we implement and define here results to be a heuristic procedure.
We conclude this section by a small illustrative example involving a network with five nodes and six edges, as shown in Fig. 2. The numbers reported on both sides of the edges indicate the passage costs through each direction for that edge. We also suppose that there are two hierarchies, for which edges (1, 2), (1, 4), (2, 5) and (4, 5) constitute hierarchy 1, whereas (2, 3) and (2, 4) constitute hierarchy 2 edges. Moreover, we also assume that the traversal costs are reduced by half after each pass. For example, the cost of the second pass through edge (1, 2) from 1 to 2 is 47 (which is the half of 94, the first passage cost).
One may notice that there are two entry nodes for hierarchy 2 which are node 2 and node 4 since these are the nodes that are incident to the edges of hierarchy 2 and hierarchy 1. Hence, graph G 0 that is constituted of the entry nodes of the network will be similar to the network given in Fig. 3. The length of edge (1, 2) in G 0 can be calculated by solving Pð1; 2; 1Þ and it is easy to see that the optimal path starting from 1 an ending at 2 that passes through each edge of hierarchy 1 is simply 1-4-5-2-1-2 with a total cost of 36 ? 30 ? 32 ? 50 ? 47 = 195. Note that the cost of passage from 1 to 2 is taken as 47 since edge (1, 2) is already traversed from 2 to 1. Similarly, the optimal path starting from 1 and ending at 4 that passes through each edge of hierarchy 1 is 1-4-5-2-1-4 with a total cost of 36 ? 30 ? 32 ? 50 ? 18 = 166. Moreover, the optimal path starting from 2 and ending at the starting node 1 that passes through each edge of hierarchy 2 is 2-3-2-4-1 with a total cost 25 ? 21 ? 33 ? 22.5 = 101.5. Note that the passage cost from node 3 to 2 is taken as 21 (half of 42, the first passage cost) since edge (2, 3) is already traversed from 2 to 3 and the passage cost from node 4 to 1 is taken also as 22.5 (half of 45) since it is already traversed during the first leg in which all the hierarchy 1 edges are passed. Besides, the optimal path starting from 4 and ending at the starting node 1 that passes through each edge of hierarchy 2 is 4-2-3-2-1 with a total cost of 18 ? 25 ? 21 ? 25 = 89. Note that passage costs from 3 to 2 and from 2 to 1 are also taken as half of the original costs since the related edges are passed before. Finally, the shortest path from i 1 ¼ f1g to i 3 ¼ f1g in graph G 0 is simply 1-4-1 with a total cost of 166 ? 89 = 255. The optimal path corresponds to 1-4-5-2-1-4-2-3-2-1 in the original network with a total cost of 255. Solving the hierarchical windy postman problem with variable service costs using a math-heuristic… 8797

Computational experiments
In this section, the selection of the parameters of the V-HWPPF and the generation of the test instances are given in the first part, and then the efficiency and accuracy of our heuristic approach are illustrated on extensive set of test instances.

Selection of the parameters and instance generation
Given that problem V-HWPP and its formulation are being introduced for the first time in this study, there are no test instances related to such problem in the existing scientific literature. We then randomly assign a number of edges to each level of hierarchy so that they sum up to the total edges number. For instance, if the number of edges is 13 and the number of hierarchies is 3, then the number of edges in hierarchies 1, 2 and 3 can take the values of 3, 6 and 4, respectively. In order to literally create the edges, we randomly select two nodes and form an edge and assign its hierarchy (starting from the first hierarchy) if the edge is connected to the previously generated edges belonging to the same hierarchy. If the created edge is not connected to the previously generated edges belonging to the same hierarchy, then we delete the edge and repeat the procedure until the number of edges belonging to the hierarchy reaches the assigned number of edges for that hierarchy. We then proceed to generate the edges of the next hierarchy. For instance, if the number of edges in hierarchy 1 is 3, then the first edge is created randomly between two randomly selected nodes and its hierarchy is assigned as 1. We then keep randomly generating the next edges until we find an edge that is connected to the first edge and assign its hierarchy as 1 too. We repeat the same procedure for edge 3 and assign its hierarchy as 1 and then proceed for generation of edges of hierarchy 2, and so on. Moreover, the first edge of a given hierarchy is forced to be connected to at least one of the higher hierarchy edges. In this way, we ensured that hierarchies are interconnected and that the edges within each hierarchy show the linear relationship. The number of hierarchies are selected as two, three, four and five, implying that we have four different number of hierarchies. The first passage (or service) cost of each edge ði; jÞ is generated randomly according to the expression c ij1 ¼ 30 þ 70 Â rand 0; 1 ð Þ where randð0; 1Þ is a uniform random number generated from ð0; 1Þ interval. c ji1 values for each edge ði; jÞ are also calculated in a similar manner. Hence, the first passage costs are random values within the interval ð30; 100Þ for all edges. The values of the other cost parameters, i.e., c ijk and c jik values for each i; j ð Þ 2 E; k [ 1, are generated using the formulas c ijk ¼ c ijðkÀ1Þ 2 þ randð c ijðkÀ1Þ 2 Þ and c jik ¼ c jiðkÀ1Þ 2 þ randð c jiðkÀ1Þ 2 Þ where randðmÞ denotes a real number randomly generated within the interval ð0; mÞ. Note that the resulting c ijk and c jik values naturally have a non-increasing structure, and they are all positive due to the employed generation mechanism.

Accuracy and efficiency of the WLA heuristic
In this section, we assess the performance of the layer algorithm by comparing the minimum distances found by the WLA heuristic with the distances and the minimum possible lower bound values reported by the state-of-the-art MILP commercial solver Gurobi (2020) on the generated test sets. We code the WLA in Visual Studio environment with C# language and we carry out all experiments using a single Intel i7-8750H core. We let Gurobi run for at most three hours for each of the test instances and the minimum objective function values and the lower bound values found in the allowed computation times are reported. However, although Gurobi is able to produce feasible solutions for all instances with 10 nodes, it cannot produce feasible solution for any instance having more than 10 nodes. In addition, Gurobi is not able to produce even a lower bound larger than zero for four instances with 10 nodes and for none of the instances having 20 nodes. On the other hand, WLA produces feasible solutions for all the instances even for the largest instances having 100 nodes. The objective function values and the lower bound values found by Gurobi, and the objective function values found by WLA are reported in Table 3 (in the Appendix). For all the instances for which Gurobi was not able to produce an upper bound and/ or a feasible solution, we used the notation NA (for nonapplicable) in the related cell.
The first and second columns of Table 3 specify the number of nodes and number of edges, respectively, while the third column reports the number of hierarchies. The fourth and fifth columns report the lower bounds and minimum cost values obtained by Gurobi while the sixth column reports the minimum cost values found by WLA. Seventh and eighth columns represent respective percentage differences between the minimum costs found by WLA and Gurobi, and between the minimum cost found by WLA and the lower bound reported by Gurobi. The percentage difference values are calculated as 100 Â c G Àc WLA c WLA and 100 Â c WLA ÀLB c WLA where c WLA and c G represent the cost values reported by WLA and Gurobi, respectively, while LB stands for the lower bound value produced by Gurobi. Finally, the nineth and tenth columns represent the computation times (in seconds) needed by Gurobi and WLA to reach the obtained solution, respectively. One can extract from Table 3 that the mathematical model results found by Gurobi outperforms WLA (by 2.24%, %0.68, and %0.47, respectively) only for three instances having 10 nodes, which are, respectively, the instances with 13 edges and five hierarchies, 18 edges and four hierarchies, and 18 edges and five hierarchies. For the one with 13 edges and five hierarchies, Gurobi's result is also the optimal solution. Both methods are able to find the optimal solutions for two instances with 10 number of nodes, and for another two instances both methods produced the same result that is higher than the lower bound found by Gurobi. Finally, WLA is better than the mathematical model (respectively, by 4.42, 18.18, 10.26, 18.63, and 13.47%) for five of the instances with 10 nodes. In summary, the objective function values of WLA and Gurobi results for the instances with 10 nodes are similar but slightly in favor of WLA. The percentage difference between the average results for the instances with 10 nodes is 5.13% in favor of WLA. More interestingly, the average computation time used by Gurobi is 8390.97 s to reach these results while WLA spends only 1.43 s on average. Hence, for the smallest instances with 10 nodes, WLA is able to produce results that are at least as successful as the model's results found by Gurobi but using much less computation time. On the other hand, Gurobi is not able to produce a feasible solution for any instance having at least 20 nodes, although it uses the entire allotted time which is three hours. We denote in the sequel as the Gurobi instances those instances for which Gurobi is able to produce a feasible solution. For a better general view of the results, we report in Table 2 the average values (extracted from Table 3) and also visualize the same objective function values for Gurobi instances in Fig. 4.
Besides, we also report the LB values found by Gurobi in order to better assess the quality of WLA results. Nevertheless, Gurobi is able to produce a lower bound larger than zero for eight instances with 10 nodes and it is not able to produce any lower bound larger than zero for the instances with at least 20 nodes. This implies that Gurobi is not able to solve even the linear relaxation of the problem for those instances within the given three hours time limit. The instances for which Gurobi is able to produce a lower bound larger than zero is called as Gurobi lower bound instances. We illustrate WLA results and the lower bound values for the Gurobi lower bound instances in Fig. 5 below.
As can be seen from Fig. 5, WLA produces almost the same LB values for the relatively small instances having 13 edges. Indeed, the average percentage difference between WLA results and the LB values are only 4.76% for the instances with 10 nodes and 13 edges. However, the average difference increases to 81.12% for the instances with 18 edges and for the remaining instances Gurobi is run (that is the instances with 10 nodes and 30 edges and the instances with 20 nodes) Gurobi reports zero as a lower bound implying that the percentage difference to be 100%. However, this does not directly indicate the deterioration of WLA performance with the increase in the instance size since even the lower bound values found by Gurobi worsen as the problem sizes increases. As a matter of fact, WLA is able to produce feasible solutions for much larger instances, reaching 100-node networks, whereas lower bounds larger than zero cannot be found by Gurobi for instances with 10 nodes and 30 edges and for all instances with 20 nodes. This may indicate that the capability of Gurobi in generating lower bounds is very sensitive to the increase in the instances size. For instance, it is possible to observe a rapid decrease in the lower bounds provided by Gurobi when the number of edges is increased from 13 to 18 for the instances having 10 nodes, whereas the objective function values of WLA possess a smooth linear-like structure increasing nearly at a constant rate. Hence, the Solving the hierarchical windy postman problem with variable service costs using a math-heuristic… 8799 increase in the percent deviation between WLA results and Gurobi's lower bounds cannot be attributed to the bad quality of WLA results but rather to the rapid deterioration in the reported Gurobi lower bound values as the instances size gets larger.
In addition, Fig. 6 reports the average computational times employed by Gurobi and WLA. One may note that the average time spent by WLA is so small for all instances with 10, 20 and 30 nodes so that they are almost at negligible level compared to the average times spent by Gurobi (this is the reason for which there is no visible average computation time of WLA for those instances). Another observation is that WLA also uses negligible amount of time (around 0.5 min) for all instances with 30 nodes on the average while it, respectively, uses around 2, 5 and 20, minutes for instances with 40, 50 and 75 nodes on the average. Finally, the average time spent for 100-node instances is around 56 min. On the other hand, the average Hence, Gurobi uses all the allotted computation time for all the instances with at least 20 nodes. In summary, it is possible to claim that not only WLA is more accurate than the mathematical model solved by the commercial solver Gurobi, in the sense that the costs of the routes generated by WLA are much lower than Gurobi instances, but also that WLA is more efficient in generating the routes in much less computation times. Finally, we observe that the running time of WLA exceeds three hours for only one instance (with 100-node instances and 1980 edges). For this reason, we avoided to test our WLA on larger instances since we think that we have achieved already significant results while solving instances that can be met in several real-life applications. The objective function values found by WLA for all the instances are depicted, for convenience, in Fig. 7, where it is clear that the objective function values increase in a smooth manner, which is another indicator of persistent achievement of WLA.

Conclusions and future research
This study dealt with the Hierarchical Windy Postman Problem with variable cost, which has not been discussed in the existing scientific literature. First, we propose a mathematical model for the problem. Our results confirmed that the exact solution procedure was not able to solve large-scale problems. For this reason, a novel heuristic method, inspired from the layer algorithm of Dror et al. (1987), has been proposed to solve the problem. Several test instances have been randomly generated in order to analyze the performance of our heuristic approach and to compare it with the performance of the mathematical model solved by the commercial solver Gurobi. The obtained results have shown that our windy layer algorithm generates high quality solutions in a very limited amount of computational time compared to the mathematical model results solved by Gurobi for the smallest sized problems with 10-nodes. For the larger sized instances, the model solved by Gurobi fails to generate feasible solutions and upper bound values, if available, are huge. In the case of 10-node instances, both the methods produce, on average, very close objective function values. However, when the average solution times are compared in the same problems group, it can be seen that the windy layer algorithm generates solution in much smaller computation times. When the size of the problems grows, Gurobi fails to produce any feasible solution and not even lower bound values for slightly larger instances whereas our heuristic method succeeds to solve instances with up to 100 nodes, 3300 edges and five hierarchies.
This work can be expanded in different ways. First, the class of problems characterized by hierarchy relations that are not linear can be discussed. Another direction of research can be analyzing the situation where service costs are fuzzy or stochastic. A variant of the problem for which the edges have different number of service requirements or whose service is characterized by a periodic nature (see Triki et al. 2017) can also be studied. Moreover, the hierarchic rural postman variant of the problem, in which only a subset of edges is required to be served, is an exciting problem that should attract the attention of researchers in future. One may also try to add some valid inequalities to formulation (20)-(32) in order to produce a stronger formulation for the open WPP aiming to have shorter computation times. Similarly, employing a heuristic to solve the open WPP would be another interesting avenue of investigation. Finally, the variants of the problem that involves multiple vehicles and/or multiple depots can also be explored.   Solving the hierarchical windy postman problem with variable service costs using a math-heuristic… 8803 Author's contribution All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by all authors. The first draft of the manuscript was written by MEK and both the other authors revised and commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Funding Open access funding provided by Università del Salento within the CRUI-CARE Agreement. None.
Data availability Enquiries about data availability should be directed to the authors.

Declarations
Conflict of interest The authors of this research certify that there is no any affiliation with or involvement in any organization or entity with financial interest or non-financial interest in the subject matter or materials discussed in this manuscript.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.