Routing for unmanned aerial vehicles: Touring dimensional sets

ARTICLE IN PRESS JID: EOR [m5G; July 26, 2021;4:30 ] European Journal of Operational Research xxx (xxxx) xxx Contents lists available at ScienceDirect European Journal of Operational Research journal homepage: www.elsevier.com/locate/ejor Discrete Optimization Routing for unmanned aerial vehicles: Touring dimensional sets Justo Puerto 1 , Carlos Valverde 1 , ∗ Department of Statistics and Operations Research, University of Seville, Seville 41012, Spain a r t i c l e i n f o Article history: Received 18 January 2021 Accepted 29 June 2021 Available online xxx Keywords: Routing Networks Logistics Conic programming and interior point methods a b s t r a c t In this paper we deal with an extension of the crossing postman problem to design routes that have to visit different shapes of dimensional elements rather than edges. This problem models the design of routes of drones or other vehicles that must visit a number of geographical elements to deliver some good or service and then move directly to the next using straight line displacements. We present two families of mathematical programming formulations. The first one is time-dependent and captures a number of characteristics of real applications at the price of using three indexes variables. The second family of formulations is not time-dependent, instead it uses connectivity properties to ensure the proper definition of routes. We compare them on a testbed of instances with different shapes of elements: second order cone (SOC) representable and polyhedral neighborhoods and polygonal chains. The computational results reported in this paper show that our models are useful and our formulations can solve to optimality medium size instances of sizes similar to other combinatorial problems including neighborhoods that have already been studied in the literature. To address larger instances we also present a heuristic algorithm that runs in two phases: clustering and Variable Neighborhood Search. This algorithm performs very well since it provides promising feasible solutions and, in addition, it can be used to initialize the solvers with feasible solutions. ©2021 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) 1. Introduction Drones, or UAVs (unmanned aerial vehicles), provide new opportunities for improving logistics in a variety of settings. Specifically, we would like to emphasize, among other characteristics, their capability for moving without an underlying network using straight line displacements. Recent technological improvements as battery life, better communication devices and reduction in manufacturing costs have increased the use of drones in logistics. Thus, this technology has increased its use in many different fields as disaster management in remote regions (see Knight, 2016 ), parcel delivery as shown in Lavars (2015) , communication coverage, worked in Amorosi, Chiaraviglio, D’Andreagiovanni, & BlefariMelazzi (2018) , traffic monitoring, infrastructure inspection, coastal surveying and many other applications. The reader is referred to the review by Otto, Agatz, Campbell, Golden, & Pesch (2018) for further references. The availability of this new technology has brought new business opportunities and, at the same time, has opened a lot of new ∗Corresponding author. E-mail addresses: [email protected] (J. Puerto), cvalv[email protected] (C. Valverde). 1 Both the authors contributed equally to this work. challenges in the Operations Research field to propose solutions to new emerging problems in the areas of logistics and routing. As drones play a growing role in business operations, questions of planning and optimization increase in practical and academic importance. However, some of the characteristics of drone’s displacement are not fully exploited by most previous routing models in literature. Unlike standard ground vehicles that must follow paths, drones can use direct connections by straight lines between destinations because they can fly across areas, but their limited battery autonomy range makes the problem of coordination with mothership vehicles a challenging problem. In 1962, Meigu Guan introduced the undirected Chinese Postman Problem (CPP) whose aim is to determine a least-cost closed route that traverses all edges of the graph. Orloff (1974) extended the CPP to travel through a subset of required edges that is known as the Rural Postman Problem (RPP). Based on this idea, Garfinkel & Webb (1999) introduced the Crossing Postman Problem (XPP) which relaxes the RPP to the case in which it is permitted to leave the edges of the network and cross from one edge to another at points other than the original vertices. These Arc Routing Problems (ARP) are studied in depth in Corberán & Laporte (2015) . On the other hand, some drone routing problems inherit some of the structure of the well-known Traveling Salesman Problem with https://doi.org/10.1016/j.ejor.2021.06.061 0377-2217/© 2021 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) Please cite this article as: J. Puerto and C. Valverde, Routing for unmanned aerial vehicles: Touring dimensional sets, European Journal of Operational Research, https://doi.org/10.1016/j.ejor.2021.0 6.0 61 J. Puerto and C. Valverde European Journal of Operational Research xxx (xxxx) xxx ARTICLE IN PRESS JID: EOR [m5G; July 26, 2021;4:30 ] Fig. 1. An example of convex sets and polygonal chains considered in the problem. neigborhoods (TSPN) that was first introduced by Arkin & Hassin (1994) and later was studied, among others, by Gentilini, Margot, & Shimada (2013) using convex sets and Yuan & Zhang (2017) presenting a hybrid framework in which metaheuristics and classical TSP solvers are combined strategically to produce high quality solutions for TSPN with arbitrary neighborhoods. Some other combinatorial optimization problems analyzed with neighborhoods are shortest paths in Disser, Mihalák, Montanar, & Widmayer (2014) , minimum spanning trees in Yang, Lin, Xu, & Xie (2007) , Blanco, Fernández, & Puerto (2017) , ordered p-median location, in Blanco (2019) and hub location, Blanco & Puerto (2021) . The aim of this paper is motivated by the design of drones’ routes that must connect a number of dimensional targets with given shapes, that we will call from now on elements , that are located on an area. The use of this terminology is not new and the interested reader is referred to Schöbel (2015) , Díaz-Báñez, Mesa, & Schöbel (2004) and Mallozzi, Puerto, & Rodríguez-Madrena (2019) for further details and references on the concept of dimensional facilities. In addition, in some cases it will be required some extra service beyond the simple visit to an element. For instance, one has to visit a percentage of its total length (assuming that its dimension is one). In our approach we would like to exploit some new features of Mixed Integer Non-Linear Programming (MINLP) to develop formulations and solution algorithms. Obviously, we have to impose some limits to the shapes of the considered elements to achieve tractable models. As a first building block, we restrict ourselves to two main types of elements (see Fig. 1 ): convex bodies and piecewise linear chains (including segments). For the case of the convex bodies, they can represent regions that the drone must reach and where the customers are willing to pick up the orders (they can be seen as uniform probability densities). On the other hand, polygonal chains can be used to model different paths that the drone must follow to do some inspection or to avoid some barriers that can appear in a real-world scenario. The use of convex bodies in the model can be extended to more general shapes as explained in Section 2 . We assume a structure of costs trying to capture the main features of these situations. We assume that there are two types of costs: 1) the travel cost of moving between elements, and 2) the travel cost of crossing (moving on) an element. The travel costs of moving between elements may change over time: it may be cheaper to go from A to B at time 1 than at time 2. On the other hand, the cost of crossing an element may be cheaper or more expensive than moving between them: controlling the drone over polygonal chains to do some inspection may be more expensive than flying directly between targets. However, one may obtain some discount for flying over some large area (parks, lakes, natural reserves...) because the drone can do a secondary job, as reporting information in its way back to the base. This fact is represented as a weighting factor in the objective function as explained in Section 2 . A survey of these coverage path planning problems can be found in Otto et al. (2018) . The goal of the model considered in this paper is to find a minimum total cost route that visits all the elements and traverses some proportions of those with dimension one. In the rest of the paper, we will refer to this problem as the Crossing Postman Problem with Neighborhoods (XPPN). The contribution of this paper is to introduce new models for the design of routes that combine several characteristics that have not been previously analyzed simultaneously: design of routes without underlying graph structure, required targets (like in the RPP) defined on dimensional elements (as in the TSP with neighborhoods) that can be polygonal chains or other kind of more general sets and free entry and exit points over the elements. Combining these features altogether gives rise to a challenging new problem that is analyzed for the first time in this paper. The paper is structured in 8 sections. The first section is the introduction. In the second section we describe the problem and set the notation followed in the rest of the paper. Section 3 is devoted to present different valid formulations of the problem. In Section 4 we present a heuristic algorithm for solving XPPN. This heuristic has two phases: clustering and Variable Neighborhood Search (VNS). The results show that it provides good quality solution in very limited computation time. Section 5 deals with some strengthening of our formulations: pre-processing variables 2 J. Puerto and C. Valverde European Journal of Operational Research xxx (xxxx) xxx ARTICLE IN PRESS JID: EOR [m5G; July 26, 2021;4:30 ] and deriving valid inequalities to be added to the formulations. Next, in Section 6 we present a decomposition algorithm ‘a la’ Benders that can be also applied to solve the problem. We derive all the details of that decomposition and show preliminary computational results. An extensive computational experience is reported in Section 7 . There, we compare the different formulations in terms of final gaps and computing time. The paper ends with a section devoted to conclusions and extensions, where we list some interesting open lines of research connected with the problems addressed in this paper. 2. Description of the problem Let V be a set of points (vertices) embedded in R 2 . (The reader may note that extensions to the three dimensional space are possible at the price of increasing the models’ complexity.) Associated with each vertex v ∈ V, we assign an element N v that can belong to two different types: either a convex set or a polygonal chain. Later, we will show how to extend the elements to deal with union of convex sets. In the former case, let C v ⊂R 2 denote the convex set associated to v that must contain v in its interior. In the latter, let P v ⊂R 2 denote the polygonal chain assigned to v that we assume to be parameterized by its breakpoints A 1 v , . . . , A n v +1 v , where n v is the number of line segments of the polygonal chain. We denote V C = { v ∈ V : v is associated with a convex set } , V P = { v ∈ V : v is associated with a polygonal chain } . Let us denote by x i v ∈ N v : i = 1 , 2 , v ∈ V the access ( x 1 v ) and exit ( x 2 v ) points to the elements N v associated with vertices v ∈ V . A feasible solution to the XPPN problem consists of a set of pairs of access and exit points, X =  v ∈ V { x 1 v , x 2 v } , together with a tour T that the drone must traverse on the graph G = (X, E) , with edge set E = E out ∪ E in , where: E out = { (x 1 v , x 2 w ) : v  = w ∈ V } , E in = { (x 1 v , x 2 v ) : v ∈ V } . Edges in the set E out are links between different elements whereas those in E in are those that define the part of the tour that is traveled within the convex neighborhoods or the polygonal chains while the drone is doing a secondary job. Observe that all the links in E in are required and therefore they must be visited by the route. Edge lengths of an outside link (x 1 v , x 2 w ) is given by the Euclidean distance, d v w (x 1 v , x 2 w ) =  x 1 v −x 2 w  2 , between their endpoints. Edge length, d v (x 1 v , x 2 v ) , of an inner link (x 1 v , x 2 v ) is computed as the distance measured over the corresponding element (polygonal or convex set). Observe that in the case of a polygonal the distance is computed as the sum of the lengths of the corresponding edges or partial edges, since the drone is following the path given by the polygonal chain. The cost of a feasible solution (X, T ) is then given by the overall sum of outside edges plus the weighted sum of the inner edges: d(X, T ) =  e v w =(x 1 v ,x 2 w ) ∈T d v w (x 1 v , x 2 w ) +  e v =(x 1 v ,x 2 v ) ∈T f v d v (x 1 v , x 2 v ) , where f v is a weighting factor for traveling within the neighborhoods. This factor depends on the worth given to a possible secondary job done by the drone. We point out that in case of overlapping of two or more neighborhoods the discount factor is accounted for each one of them, as shown in the above formula. The reader may note that in all our discussions we are assuming that the autonomy of the drone battery suffices to travel the whole route. Therefore, the model does not allow a route longer than the flying autonomy. Throughout this paper we adopt the following notation: •T G as the set of incidence vectors associated with tours on G, i.e., T G = { z ∈ R | E| + : z is a tour on G } . •X =  v ∈ V ( N v ×N v ) , the space where the access and exit points are selected. The goal of XPPN is to find a feasible solution (X, T ) of minimal total cost. Then, it can be expressed as: min  e v w =(x 1 v ,x 2 w ) ∈ E out d v w (e ) z e +  e v =(x 1 v ,x 2 v ) ∈ E in f e d v (e ) s.t z ∈ T G , x ∈ X (1) Here it is assumed that the drone route enters and exits from an element only once. Note that, since the distance between neighborhoods is minimized, there always exists an optimal solution in which the drone visits each neighborhood only once. The reader may observe that the above formulation is only formal, but it is clearly not separable into the continuous and discrete counterparts since the access and exit point to each one of the elements (continuous part) depend on the order of the visit to the elements (discrete part) and vice versa. We also point out that the discrete part, that is a TSP, is an NP-hard problem whereas the continuous part, that is a location problem, is easily solvable by using interior-point algorithms. This structure is exploited to decompose the problem in a master problem (TSP) and a subproblem (Location Problem) in the Benders decomposition (see Section 6 ). Moreover, the problem involves Euclidean distances among variable points and sets, therefore it is not linearly representable. In spite of that, it is suitable to model this problem as a MINLP. In this paper, we focus on the case where the sets C v are second order cone (SOC) representable, that is, the sets can be expressed by using second-order cone constraints as follows: x i v ∈ C v ⇐⇒  A j v x i v + b j v  ≤(c j v ) T x i v + d j v , j = 1 , . . . , n v , ( C −C ) where x i v , i = 1 , 2 is the decision variable, A j v , b j v , c j v and d j v are parameters of the constraint jand n v represents the number of constraints that appear in the block associated to vertex v . Note that these inequalities can also model linear constraints (for A j v , b j v ≡0 ), ellipsoids and hyperbolic constraints (see Lobo, Vandenberghe, Boyd, & Lebret, 1998 for more details). These type of elements could be extended further to unions of SOC representable sets. This type of neighborhood is obtained introducing binary variables, whose meaning is similar to those in disjunctive programming. Thus, we can determine in which set of the union happens the access or the departure points of the different sets. Let {C 1 v , . . . , C m v v } be the second order cone representable sets that define the neighborhood associated to the vertex v and let U v = m v   =1 C  v denote the union of these sets. Consider the binary variable χi v that assumes the value of one if x i v is located in the set C  v and zero otherwise. Thus, for each v ∈ V, we can model that x i v ∈ U v by using the following inequalities for each i = 1 , 2 : x i v ∈ U v ⇐⇒   A j  v x i v + b j  v  ≤(c j  v ) T x i v + d j  v + M j  v (1 −χi v ) ,  = 1 , . . . , m v , j  = 1 , . . . , n  v ,  m v  =1 χi v = 1 , ( U −C ) where M j  v is a big-M constant on the maximal distance between two points in the union of sets. The reader may observe that one can replace ( C-C) by ( U −C ) in all our formulations without compromising their validity. Therefore, our model can deal easily with these more general forms of neighborhoods. On the other hand, the second type of elements are the piecewise linear constraints. Let n Sv be the number of line segments of the polygonal chain v . Since we need to refer to interior points of the segment, these continuum of points is parametrized by the two endpoints of the segment: x ∈ [ A j v , A j+1 v ] if and only if 3 J. Puerto and C. Valverde European Journal of Operational Research xxx (xxxx) xxx ARTICLE IN PRESS JID: EOR [m5G; July 26, 2021;4:30 ] ∃ γ∈ [0 , 1] such that x = γA j v + (1 −γ) A j+1 v . In order to deal with them, we introduce the following variables for each vertex v ∈ V P and i = 1 , 2 : •u v : Binary variable that determines the traveling direction in the polygonal chain v . •γij v : Continuous variable in [0,1] that represents the parameter value of the x i v variable in the line segment jof the polygonal chain v , j = 1 , . . . , n Sv . •μij v : Binary variable that is one when x i v is located in the line segment jof the polygonal chain v , and zero otherwise, for j = 1 , . . . , n Sv . •λi v : Continuous variable in [0 , n Sv ] that models the parametrization of the entry or exit points along the polygonal chain associated with v . Using these variables, we can determine the placement of the entry and exit points on the polygonal chain v introducing the following inequalities for each i = 1 , 2 : x i v ∈ P v ⇐⇒ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ λi v −j ≥γij v −(n Sv + 1)(1 −μij v ) , j =2 , . . . , n Sv + 1 λi v −j ≤γij v + (n Sv + 1)(1 −μij v ) , j =2 , . . . , n Sv + 1 γi 1 v ≤μi 1 v γij v ≤μij−1 v + μij v j =2 , . . . , n Sv γin Sv v ≤μin Sv v  n Sv j=1 μij v =1  n Sv +1 j=1 γij v =1 x i v =  n Sv +1 j=1 γij v A j v ( P −C ) Observe that the first and second inequalities determine the upper and lower limits for the parametrization of each segment of P v . If μij v = 0 the inequalities are always fulfilled and there is no entry or exit point in the jth segment of the polygonal v . On the contrary, if μij v = 1 then λi v ∈ [ j, j + 1] meaning that the corresponding entry or exit point is in the jth segment of the polygonal P v . The third, fourth and fifth inequalities link μij v and γij v variables (and thus implicitly λi v ):they state that the variable γij v that gives the representation of a point x i v on the line segment jis active (non-null) only if this line segment is chosen (to enter of exit), i.e., μij v = 1 . The sixth equation sets that only one line segment is chosen for entering or leaving each polygonal chain. Finally, the seventh equation and eighth inequality set the representation of x i v as a convex combination of the extreme points of the adequate line segment. In addition, we assume that the tour must traverse at least some given percentage αv of each polygonal chain total length. Denoting by λmin v and λmax v the parameter values of λrepresenting the initial and final points of P v , respectively, we can model that condition by the following absolute value constraint: | λ1 v −λ2 v | ≥n S v αv ⇐⇒ ⎧ ⎪ ⎨ ⎪ ⎩ λ1 v −λ2 v = λmax v −λmin v λmax v + λmin v ≥αv n Sv λmax v ≤n Sv (1 −u v ) λmin v ≤n Sv u v . (α−C ) The above modelling assumptions are sufficient to address the range of situations that we want to model. Obviously, they could be more general at the price of not being easy to implement with off-the-shelf solvers. 2.1. Some interesting particular cases Three very interesting well-known models appear as particular cases of the problems that can be modelled within our framework. If the element associated with each vertex is a single point the problem reduces to the standard traveling salesman problem. If the element associated with each vertex v ∈ V is a segment P v = [ x 1 v , x 2 v ] and αv = 1 , then XPPN becomes the classical Rural Postman Problem in which the edges (x 1 v , x 2 v ) are required, in the Fig. 2. An example with 9 elements: 7 convex sets and 2 polygonal chains. complete graph induced by these vertices with edge lengths given by the Euclidean norm distance (see Orloff, 1974 ). In addition, if αv  = 1 , ∀ v ∈ V then XPPN is an extension of the RPP where some edges are only partially required. On the other hand, if the considered neighborhoods are big enough so that ∩ v ∈ V C v  = ∅ , then the problem reduces to finding a degenerate one-vertex tour and the solution to the XPPN is that vertex with cost 0. Finally, if all elements N v , v ∈ V are neighborhoods we obtain the Traveling Salesman problem with Neighborhoods (see Arkin & Hassin, 1994 ). Fig. 2 shows an example of the solution obtained for a case in which the elements are circles, triangles and we also have two polygonal chains to visit in our required route. The discussion above allows us to state the complexity of the XPPN. Theorem 1. The decision version of the problem XPPN, given a length L deciding whether the graph G has a XPPN tour of length at most L, is NP-complete. The proof follows using a reduction from TSP that as shown above is a particular case of this problem. 3. Mixed integer non linear programming formulations In this section we present alternative MINLP formulations for the XPPN that will be compared computationally in Section 7 . First, we start with a time dependent formulation that allows us to include a number of specific characteristics in the modeling phase such as time dependent travel distances, time windows or time dependent discount factors. Then, we give another formulation that does not make reference to stages in the routes and that simplifies the model at the price of losing some of the above mentioned characteristics. 3.1. A time dependent formulation One way to model the drone route in our problem is to make variables dependent on the index of the stage when an element is visited in the sequence of visited elements. Thus, this formulation requires binary variables depending on the index order when 4 J. Puerto and C. Valverde European Journal of Operational Research xxx (xxxx) xxx ARTICLE IN PRESS JID: EOR [m5G; July 26, 2021;4:30 ] they are chosen. Since variables depend on time parameters in the problem, weighting factors for visiting the neighborhoods (b t v ) and distances d t (as proxy for travel times βt ), can also be dependent on the stage when they are used. To model the problem, we introduce a binary variable y t v to indicate that the element associated with vertex v is visited by the drone at stage t. In addition, we define the following variables: •y t v : Binary variable whose value is one when v is visited at the tth position in the route sequence and zero otherwise. •z t v w : Binary variable that is one when v and w are visited consecutively, assuming that v is visited at the stage tand zero otherwise. •z t v w = y t v y t+1 w , v  = w . •d t v w : Continuous variable that represents the distance between pairs of chosen points v , w from different components at the stage t. •d t v : Continuous variable that represents the distance between two consecutive points within the same component associated with v ∈ V at the stage t. •λ1 v , λ2 v : Continuous variables determining the position of x 1 v and x 2 v , respectively, in the polygonal chain P v . Using these variables, the first formulation follows: min | V |  t=1  v  = w d t v w z t v w + | V |  t=1  v ∈ V f t v d t v (2a) s.t. d t v w ≥βt v w  x 2 v −x 1 w  , ∀ v  = w (2b) d t v ≥βt v  x 1 v −x 2 v  , ∀ v ∈ V (2c)  v ∈ V y t v = 1 , ∀ t (2d) | V |  t=1 y t v = 1 , ∀ v ∈ V (2e) y t v + y t+1 w −1 ≤z t v w , ∀ v  = w, t = 1 , . . . , |C| −1 (2f) (C −C ) , (P −C ) , (α−C ) (2g) The first addend of the objective function (2a) includes the drone traveling distance among different elements while the second one accounts for the distances between the entry and exit points of each component taking into account the weighting factor for traveling within this component at the stage t. Constraints (2d) and (2e) state, respectively, that in each stage the route visits one element and each component is traversed once and only once. Constraint (2f) is obtained by linearizing z t v w and ensures that if we travel from v to w, assuming that we are in v at the instant t, then we visit v in tand w in t + 1 . Constraint (2g) refers to the domain of the entry and exit points of each element in the problem, as well as the minimal required percentage of the polygonal chain length that must be traversed by the drone. They were defined in Section 2 . Despite the versatility of this formulation for capturing actual characteristics of drone routes, its drawback comes from the three index dimension of its variables which makes it difficult to handle medium size instances. In the next section, we shall simplify this formulation making it independent of time at the price of losing some of its time-dependent characteristics. 3.2. Non-time dependent formulations The simplification mentioned above can be performed, based on the rationale of ensuring connectivity on the graph G, through different sets of inequalities. In particular, we compare MillerTucker -Zemlin (MTZ) inequalities and subtour elimination constraints (SEC). All formulations use the following sets of decision variables: •Binary variables z e ∈ { 0 , 1 } , e ∈ E out , to represent the edges of the tours. •Continuous variables d e ≥0 , e = { v , w } ∈ E out ⊆E, to represent the distance d v w (x 1 v , x 2 w ) between the pairs of selected points of different elements (neighborhoods) and d v ≥0 , v ∈ V, to represent the distance d v (x 1 v , x 2 v ) between the pairs of points of the same element. Let D e = { d ∈ R | E out | + : d e ≥d v w (x 1 v , x 2 w ) , ∀ e = (v , w ) ∈ E out , x ∈ X } , D v = { d ∈ R | E in | + : d v ≥d v = (x 1 v , x 2 v ) , ∀ v ∈ V, x ∈ X } , denote the domains for the feasibility of the dvariables. The reader can see that these sets, namely D e and D v , can be alternatively described using the constraints  x 1 v −x 2 w  2 ≤d e , ∀ e = { v , w } ∈ E out , ( D 1 )  x 1 v −x 2 v  2 ≤d v , ∀ v ∈ V, ( D 2 ) x ∈ X , ( D 3 ) which set the distance values and impose that x belongs to its suitable neighborhood. Then, a generic bilinear formulation for XPPN is min  e ∈ E out d e z e +  e ∈ E in f v d v ( Pdz ) s.t. z ∈ T G , ( D 1) , ( D 2) (C −C ) , (P −C ) , (α−C ) The reader should observe that, as already mentioned, the above formulation is bilinear since the first term of the objective function contains products of variables of the form d e z e , for e ∈ E out . Next, we use McCormick’s envelopes ( McCormick, 1976 ) for the linearization of those bilinear terms of the objective function. We define additional variables p e ≥0 , e ∈ E out that stand for that product. Replacing the products by the new variables and introducing a new set of constraints enforcing the correct representation, we obtain the following formulation: min P =  e ∈ E out p e +  v ∈ V f v d v ( RL-XPPN ) s.t. p e ≥d e −M e (1 −z e ) ∀ e ∈ E out ( LIN-Mc ) p e ≥0 , ∀ e ∈ E out z ∈ T G ( D 1) , ( D 2) (C −C ) , (P −C ) , (α−C ) Here M e denotes an upper bound of the distance between the sets that are joined by e . Furthermore, this formulation can be reinforced by adding some valid inequalities: p e ≥m e z e , ∀ e ∈ E out and d v ≤M v , ∀ v ∈ V, where m e and M v are bounds that are adjusted in Section 5 . The first family of valid inequalities sets lower bounds on the values for p e 5 J. Puerto and C. Valverde European Journal of Operational Research xxx (xxxx) xxx ARTICLE IN PRESS JID: EOR [m5G; July 26, 2021;4:30 ] whereas the second ones sets upper bounds on the distances traveled by the drone within each neighborhood. The above discussion leads us to strengthen a generic formulation for XPPN. This formulation will be particularized once the connectivity condition of the solutions is specifically introduced in the model. min P =  e ∈ E out p e +  v ∈ V f v d v s.t. p e ≥d e −M e (1 −z e ) ∀ e ∈ E out ( LIN-Mc ) p e ≥m e z e ∀ e ∈ E out ( VI-1 ) d v ≤M v ∀ v ∈ V ( VI-2 ) z ∈ T G ( D 1) , ( D 2) (C −C ) , (P −C ) , (α−C ) The two formulations that we present below differ from one another in the family of constraints used to enforce connectivity. One of them is by the family of subtour elimination constraints (SEC), Edmonds (2003) . The other one relies on a compact formulation based on the well-known Miller-Tucker-Zemlin (MTZ) constraints, Miller, Tucker, & Zemlin (1960) . 3.2.1. A valid formulation for XPPN based on SECs The family of SEC is well-known in combinatorial optimization. It enforces connectivity by imposing that the number of edges among any subset of vertices can not exceed its cardinality minus one. Augmenting these constraints into the generic formulation presented above we obtain the following valid formulation for XPPN: min P =  e ∈ E out p e +  v ∈ V f v d v ( SEC-XPPN ) ( LIN-Mc ) , ( VI-1 ) , ( VI-2 )  w ∈ V \{ v } z v w = 1 , ∀ v ∈ V ( C 1 )  w ∈ V \{ v } z w v = 1 , ∀ v ∈ V ( C 2 )  e =(v ,w ): v ,w ∈ S z e ≤| S| −1 , ∀ S  V ( SEC ) ( D 1) , ( D 2) (C −C ) , (P −C ) , (α−C ) Assignment Constraints ( C 1 ) and ( C 2 ) ensure that the drone enters and exits each component of the problem exactly once. Constraint (SEC) prevents the existence of subtours. This constraint forces that in any subset Sof nodes included in V there can not be more edges between nodes in Sthan its number of nodes minus one, thus avoiding the existence of cycles. Since there is an exponential number of SEC constraints, when we implement this formulation we need to perform a row generation procedure including constraints, whenever they are required, by a separation oracle. To find SEC inequalities, as usual, we search for disconnected components in the current solution. Among them, we choose the shortest subtour found in the solution to be added as a lazy constraint to the model. If the considered distance between components is symmetric, we obtain the symmetric formulation based on SECs, denoted by (sSEC-XPPN). In this formulation, we can halve the number of binary variables and replace constraints ( C 1 ) and ( C 2 ) in ( SEC-XPPN ) by the following connectivity restrictions:  w ∈ V \{ v } z w v = 2 , ∀ v ∈ V. 3.2.2. XPPN formulation based on the Miller-Tucker-Zemlin inequalities This section addresses an alternative formulation that results replacing SEC inequalities by the so called Miller-Tucker-Zemlin constraints (see Miller et al., 1960 ). In this formulation, we introduce the integer variable s v to generate an alternative formulation that eliminates the subtours and the exponential number of inequalities of ( SEC-XPPN ) . min P =  e ∈ E out p e +  v ∈ V f v d v ( MTZ-XPPN ) s.t. ( LIN-Mc ) , ( VI-1 ) , ( VI-2 )  w ∈ V\{ v } z v w = 1 , ∀ v ∈ V ( C 1 )  w ∈ V\{ v } z w v = 1 , ∀ v ∈ V ( C 2 ) | V | z v w + s v −s w ≤| V | −1 , ∀ e = (v , w ) ∈ E out ( MTZ 1 ) s 1 = 1 ( MTZ 2 ) 2 ≤s v ≤| V | , ∀ v ∈ V ( MTZ 3 ) s v −s w + | V | z w v ≤| V | −1 , ∀ e = (v , w ) ∈ E out , w > 1 ( MTZ 4 ) s v −s w + (| V | −2) z w v ≤| V | −1 , ∀ e = (v , w ) ∈ E out , v > 1 ( MTZ 5 ) ( D 1) , ( D 2) (C −C ) , (P −C ) , (α−C ) Again constraints ( C 1 ) and ( C 2 ) require that in each feasible solution only one edge departs from node v and only one edge enters at node v for any v ∈ V, respectively. It is well-known that constraints ( MTZ 1 ) - ( MTZ 3 ) (see Miller et al., 1960 ) model the elimination of subtours. The constraints ( MTZ 1 ) - ( MTZ 3 ) enforce connectivity, i.e., that there is only a single tour covering all vertices. The constraints ( MTZ 4 ) and ( MTZ 5 ) define the intermediate conditions for the tour that may improve the performance of this formulation over the formulation based on subtour elimination constraints (see Sawik (2016) for more details). Now we state a result related to the relationship between the SEC and MTZ polytopes of our formulations of the XPPN, that is, the feasible regions of the respective LP relaxations of these models. Theorem 2. The SEC polytope is contained in the MTZ polytope for the XPPN. Proof. Observe that the only difference between these two polytopes is the family of constraints that ensures the elimination of subtours. Therefore, it is enough to see that the (SEC) constraints are stronger than those given in ( MTZ 1 ) - ( MTZ 3 ) which is proved in Velednitsky (2017) .  4. A heuristic algorithm for XPPN In this section we present a heuristic algorithm for solving XPPN. This algorithm has two different applications. On the one hand, it provides good quality feasible solutions for XPPN that become a promising alternative to exact methods whenever the size of the problems is large. On the other hand, it also helps in solving exactly XPPN by feeding the exact formulations with a good initial solution which in turns speeds up the branch and bound search. The considered algorithm is composed by two phases: the Clustering Phase and the Variable Neighborhood Search (VNS) Phase. The so called clustering phase determines some points in each dimensional element (polygonal chain or neighborhood) and then the VNS phase finds a heuristic tour on the complete graph spanned by the previously obtained points. The clustering phase The first phase of the heuristic algorithm is based on solving a relatively easy single facility location problem: the Weber or median Problem. The solution of this problem looks for a prototype point (a representative) x v , v ∈ V of the dimensional elements in the problem (neighborhoods and polygonal chains) and another 6 J. Puerto and C. Valverde European Journal of Operational Research xxx (xxxx) xxx ARTICLE IN PRESS JID: EOR [m5G; July 26, 2021;4:30 ] Fig. 3. Illustration of the first phase of the heuristic algorithm. point, Med, so that the sum of the Euclidean distances from x v to Med is minimized. min  v ∈ V  x v −Med  ( Web er ) s.t. (C −C ) , (P −C ) , (α−C ) The idea of this approach is to find some points that are likely to be close to the true chosen points in each element in the final optimal drone route. Fig. 3 shows an example that combines six neighborhoods and four polygonal chains. Red points represent the points of each set and the green point is the proposed 1-median obtained after solving the corresponding Weber problem described above. The variable neighborhood search phase Once the points of each set have been chosen, the idea is to find the minimal cost drone route that joins these points. To obtain this route, we have used the well-known and general Variable Neighborhood Search metaheuristic developed in Mladenovi ´ c & Hansen (1997) . The Python implementation code has been taken from Pereira (2018) . In that implementation, the distance matrix is computed by taking the Euclidean distances between each pair of points. In our case, we had to modify it because our distance matrix requires also distances computed along the different considered polygonals. Using the example depicted in Fig. 3 , we generate a tour considering this VNS approach with a maximum number of 25 attempts, a neighborhood of size 5 and 10 iterations. The final result is shown in Fig. 4 . Finally, in order to build a feasible solution for XPPN we take into account the position of the points (represented by x 1 v and x 2 v ) and the order in which they are visited in the tour obtained by the VNS phase of our heuristic (represented by z v w ). Once the solution is built, it can also be taken as an initial solution for any of the exact formulations presented above. In the following, we present the pseudo-code of this heuristic: Fig. 4. Application of the VNS phase to the example of Fig. 3 . Algorithm 1: Heuristic for solving XPPN. Let {N v : v ∈ V } be the neighborhood set.Set at tempt s = 25 , neigh _ size = 5 , iter = 10 . 1. Solve the Weber problem for N v to get ¯ x . 2. Consider the VNS approach with parameters at tempt s, neigh _ size and iter and points ¯ x to obtain the order of visit to the neighborhoods ¯ z . 5. Strengthening the formulation of XPPN 5.1. Pre-processing In this section we explore the geometry of the neighborhoods that appear in the problem to fix a priori some variables and to increase the efficiency of the model. First of all, we consider two special cases that relate the position of the entry and exit points of each neighborhood with the coefficient f v of the objective function. Remark. If the problem verifies that f v = 0 for all v ∈ V C , then the entry and exit points x 1 v and x 2 v selected in each neighborhood are the same that the ones obtained by minimizing the distance between the neighborhoods. Remark. If f v ≥1 for some v ∈ V C , then, there exists an optimal solution verifying x 1 v = x 2 v . Proof. Let us consider an optimal route and let p be the path in that route that visits C v . Assume without loss of generality that to visit C v , the route departs of the previous element C u from x 2 u , enters C v through x 1 v and exits from x 2 v where x 1 v  = x 2 v and f v ≥1 . Assume again without loss of generality that after visiting C v the route goes to C w entering by x 1 w . Let us consider the alternative path p  formed by x 2 u , the midpoint x  v between x 1 v and x 2 v and x 1 w . The contribution of visiting C v in the objective function of the problem will be 7 J. Puerto and C. Valverde European Journal of Operational Research xxx (xxxx) xxx ARTICLE IN PRESS JID: EOR [m5G; July 26, 2021;4:30 ] length (p  ) = d(x 2 u , x  v ) + d(x  v , x 1 w ) ≤d(x 2 u , x 1 v ) + d(x 1 v , x  v ) + d(x  v , x 2 v ) + d(x 2 v , x 1 w ) = d(x 2 u , x 1 v ) + d(x 1 v , x 2 v ) + d(x 2 v , x 1 w ) ≤d(x 2 u , x 1 v ) + f v d(x 1 v , x 2 v ) + d(x 2 v , x 1 w ) = length (p) , but pis an optimal path to visit those elements in this optimal solution which turns all the above inequalities into equalities. Therefore, the path p  is also an optimal path within that optimal solution. However, by construction on C v , p  has the same entry and exit point which proves the claim.  From now on, we assume in the rest of this section that f v ≥1 for all v ∈ V . The following outcome restricts the domain where the selected points can be located. Proposition 1. There exists always an optimal solution of the XPPN whose selected points are placed in the boundary of the neighborhoods. Proof. If the number of elements of the problem is two, the problem consists of calculating the minimum distance between two convex sets and it is known that the selected points are clearly located in the boundary of the sets if they do not overlap and can be chosen in the border in case of overlapping. If the number of neighborhoods is more than two, we can reduce the proof to analyze three consecutive elements. Let T be the triangle spanned by the point x u ∈ C u of the previously visited neighborhood, the point x v in the neighborhood C v and the point x w ∈ C w of the next neighborhood to be visited in an optimal sequence. We could have three possible cases depending on the number of points allocated in the boundary. If x u , x v , x w are aligned, for each point that is not in the boundary, namely C v , we can consider the closest point obtained by the intersection of the line generated by x u and x v with the boundary of C v , ∂C v . This point is also aligned with the others and its contribution to the objective function is the same as the one given by x v . Therefore, let assume that these points are not aligned. We have three cases: •Case 1 : Suppose that only x v ∈ C v is not in the boundary. Let x u x w be the line segment that joins x u and x w . Let suppose, for the sake of contradiction, that there exists a neighborhood C v whose selected point in the optimal sequence is in the topological interior of C v , i.e., x v ∈ int(C v ) . The idea of the proof is to find another point in the boundary of C v closer to x u and x w . We consider the point x  v = r ⊥ ∩ ∂C v , where r ⊥ is the perpendicular line to the line segment x u x w and ∂C v is the boundary of C v . Observe that this intersection produces two points in ∂C v : among them we take the closest one to x u x w . If we call T  the triangle generated by x u , x  v and x w , then the height of T  to x u x w is smaller than the one of T . Hence by the Pythagoras Theorem, x u x  v < x u x v and x  v x w < x v x w , which is a contradiction. •Case 2 : Assume that x u ∈ C u and x v ∈ C v are not in the boundary. We can take x  u = x u x v ∩ ∂C u . This point is closer to x v than x u and it is in the boundary of C u . Therefore, we have two points in the boundary and we can apply the previous case to conclude that x v must be in ∂C v too. •Case 3 : Finally, suppose that no point is in the boundary of each neighborhood. Again, we can construct x  u = x u x v ∩ ∂C u that is closer to x v and x w . Then, we have a point in the boundary and Case 2 can be applied to the rest of points.  The special case in which all the neighborhoods are circles, allows us to limit even more the location of the points based on the construction given in the Proposition 1 . Corollary 2.1. Any point selected in an optimal solution of the XPPN when all the neighborhoods are circles is placed in some arc of one of the circumferences inside of the convex hull generated by the center of the circles. Proof. If we have two neighborhoods, the selected points are located in the line segment that joins the center of the circles and the result follows. If the number of neighborhoods is more than two, we can reduce the proof to analyze three consecutive elements. Let T be the triangle spanned by the point x u ∈ C u of the previously visited neighborhood, the point x v in the neighborhood C v and the point x w ∈ C w of the next neighborhood in an optimal sequence. Let assume that we have two points inside the convex hull and x v ∈ C v does not satisfy this property. Let also x u x w be the line segment that joins x u and x w . We distinguish two cases depending on the location of x v in the neighborhood: •If x u , x v , x w are aligned, it is straightforward to conclude that x v is in the convex hull of the centers. •If x u , x v , x w are not aligned, let assume that in N v its selected point is not in the convex hull C v of the centers of the neighborhoods. The idea of the proof is to find another point in the boundary of the convex hull C v whose distance to x u and x w is smaller than the distance from x v . We split the boundary of C v (circumference) in two arcs A R and A B . These arcs are built by taking the perpendicular line to the edge of the convex hull C: •If x v ∈ A R , we take x  v the projection to the convex hull and it produces a triangle T  with lower height to x u x w . Then, we can use the Proposition 1 to construct a point in the boundary of C v that lies in the convex hull. (See Fig. 5.1). •If x v ∈ A B , we construct x  v the diametrically opposite point of x v in N v . This point also produces a smaller height that contradicts the assumption that x v gives the shortest tour. (See Fig. 5.2) If the number of points outside the convex hull is more than one, we can apply this procedure iteratively to include these points in the convex hull generated by the center of the circles.  8 J. Puerto and C. Valverde European Journal of Operational Research xxx (xxxx) xxx ARTICLE IN PRESS JID: EOR [m5G; July 26, 2021;4:30 ] Finally, we conclude this section giving another result that allows one to eliminate some neighborhoods and thus simplify the problem without modifying the objective value of the problem. Proposition 2. Given two neighborhoods A and B, if B ⊃A, then B can be removed in the problem. Proof. Starting from the optimal solution of problem without B, we are going to build an optimal solution including B that is essentially the same. Let z ∗the optimal tour by deleting the neighborhood B in the problem. By connectivity, there exist two neighborhoods A −1 and A +1 that are connected with A, i.e., such that z ∗ A −1 A = z ∗ AA +1 = 1 . In addition, let x ∗ A be the point chosen to visit the neighborhood A . If we include B ⊃A in the problem and we fix x B = x ∗ A and z AB = z BA +1 = 1 . This solution is also a simple path whose objective value is the same because d(A, B ) = 0 and d(B, A +1 ) = d ∗(A, A +1 ) .  5.2. Valid inequalities The different models that we have proposed include in one way or another big-M constants. In order to strengthen the formulations we provide good upper bounds for those constants. In this section we present some results that adjust them for each kind of set considered in our models. The first big-M constant we need to adjust is M e that denotes an upper bound of the distance between the sets joined by an edge e ∈ E out . We have three cases that depend on the shape of the sets A and B : •If A and B are both ellipsoids, we cannot easily compute the maximum distance between A and B, but we can generate an upper bound of this distance by taking diametrically opposite points of minimum radius circles containing each ellipsoid. When both ellipsoids are circles, this bound coincides with the maximum distance. •If A is an ellipsoid and B is a polygon or a polygonal chain, we can set this bound by the maximum of the distances of each vertex of B to the center of A plus the radius of the minimum circle that contains the ellipsoid A . •If A and B are both polygons or polygonal chains, this bound can be computed exactly by taking the maximum of the distances between vertices of A and B . The second bound to be adjusted is m e . It denotes a lower bound of the distance of the sets joined by the edge e ∈ E out . In this case, we can compute this distance exactly by solving a convex program that minimizes the distance between the sets A and B . In the Figs. 5 –7 we show the selected maximal (red) and minimal (blue) bounds depending on the shape of the sets. In addition, the third bound represents the maximal distance between two points within a given neighborhood. We can compute this upper bound according to the shape of this set (see Fig. 8 ): •If the set is an ellipsoid, we can take diametrically opposite points of the minimum radius circle that contains this ellipsoid. Fig. 5. Upper and lower bound when both sets are ellipsoids. Fig. 6. Upper and lower bound when a set is a polygon and the other is an ellipsoid. Fig. 7. Upper and lower bound when both sets are polygons. Fig. 8. Upper bound on the maximal distance within a set. •If the set is a polygon, we can compute the maximum of the distances between each pair of vertices. •If the set is a polygonal chain, this bound equals the length of the polygonal. 6. A decomposition algorithm In this section we present an alternative row generation approach to solve the XPPN based on a Benders decomposition of the problem. The general method is based on the following observation: If we fix z ∈ T G in the generic formulation of XPPN, we obtain a continuous SOC problem, which is well-known to be convex. On the other hand, the objective function that we are considering is bilinear. Hence, we can use a Benders-like decomposition approach (see Benders, 1962 ) to generate an iterative algorithm that solves this problem. For a given ¯ z ∈ T G , the “optimal” vertices and distances of its associated XPPN can be computed by solving the following sub9 J. Puerto and C. Valverde European Journal of Operational Research xxx (xxxx) xxx ARTICLE IN PRESS JID: EOR [m5G; July 26, 2021;4:30 ] Fig. 12. Final gap after 7200 seconds. a comparison between formulation MTZ with and without the initial solution provided by the VNS heuristic. The results are summarized in Table 1 . This table reports average results for instances of sizes 5,10 and 15 neighborhoods with all combinations of radii and modes. It compares the final gap and running times for the formulation with and without initial solution ( Final Gap (Init), Opt. Time (Init) ) (resp. Final Gap (NoInit), Opt. Time (NoInit) ). The results are also depicted in the boxplox diagrams in Fig. 10 . Both, table and figure, clearly show that loading an initial solution helps in reducing the gap and the cpu time: all the instances up to 10 neighborhoods are solved to optimality with and without initial solution but for 15 neighborhoods the final gap in the first case is always better than in the latter (blue boxes are always below orange ones). Based on this results, in the following, we have always solved the instances loading an initial solution. 7.5. Comparing benders cuts with the MTZ formulation Here, the decomposition algorithm described in Algorithm 2 is compared with the MTZ formulation without initialization. The computational results obtained by our implementation are included in Table 2 . This table compares the results of the Final Gap, cpu time and number of cuts added applying the decomposition algorithm versus those obtained with formulation MTZ. From these results we conclude that the decomposition algorithm performs worse than formulation MTZ even for small size problems. The Benders optimality cuts involve big-Ms, which in turns implies that a lot of cuts are needed (if not all) to certify optimality. The big-M constraints comes from the linearization of bilinear terms which do to allow to apply the Benders approach because the lack of convexity. Thus, this may be one of the reasons why its performance is worse than MTZ. To reinforce our observation, we have also included a performance profile of number of solved instances versus time for formulations sSEC, MTZ and the decomposition algorithm (see Fig. 11 ). The reader can observe that the number of solved instances within the time limit is approximately one half comparing Benders decomposition with MTZ and sSEC. These results lead us to not include this algorithm in the final computational experience for larger problem sizes presented in the last subsection. 16 J. Puerto and C. Valverde European Journal of Operational Research xxx (xxxx) xxx ARTICLE IN PRESS JID: EOR [m5G; July 26, 2021;4:30 ] Fig. 13. Performance profile: Time vs #Solved. 7.6. Comparing MTZ, SEC and sSEC with initialization from a heuristic solution The remaining information of our computational experiments can be found in Tables 3–5 . The first one reports our results for formulation SEC, the second one for sSEC (symmetric version of SEC) and the third one for MTZ. Information in all the three tables is organized in the same way. Each row shows averages of five instances for different combinations of factors ( Size , Radii and Mode ) each one with four different levels. Our tables have 9 columns. The first three (Size, Radii and Mode) describe the parameters of the problem. Then, we report the final gap ( % Final Gap ), time required by the exact method ( Exact Time ), time required by the heuristic ( Heur. Time ) and the % improvement of the gap with respect to the initial solution ( % Improved Gap ). Overall, we have solved 320 instances. To have a clearer view of the results we also present some comparative boxplots obtained from the tables above. First of all, we report the final gap after two hours of running time. We have gathered all the information in Fig. 12 . It is organized in four rows corresponding with the different modes: row i shows results for Mode i, i = 1 , . . . , 4 . Within each row, there are four columns one per radius size. Then, each graph within this 4 ×4 grid contains comparative diagrams for the four different problem sizes considered, namely | V | = 5 , 10 , 15 , 20 neighborhoods. Finally, for each problem size we compare the results obtained for our three different formulations Miller-Tucker-Zemlin (MTZ), Subtour elimination (SEC) and the symmetric version of SEC (sSEC). For instance, looking at the second row, third column (Mode 2, Radius 3) one can see that for | V | = 5 and 10, which correspond to the first two boxes the gap of the three formulations is zero in all the instances (actually the boxes are collapsed to lines). However, for | V | = 15 and 20 MTZ seems to outperforms SEC and sSEC, and moreover, sSEC is also better than SEC since the green boxes lie below the orange ones. As a general comment, one can observe that for all combinations of factors MTZ (the blue boxes) outperforms SEC (orange) and sSEC (green) and also sSEC reports smaller gaps than SEC, with the only exception of Mode 3 where SEC seems to work better than sSEC. Finally, we have included in Fig. 13 a performance profile graph of number of instances solved versus time. This figure shows that SEC formulation is the weakest since it solves less number of instances in the same cpu time. The comparison between MTZ and sSEC is not that clear although in the long run MTZ seems to outperforms sSEC since the former solves more instances than the latter. We also compare next, the behaviour of SEC and sSEC in number of cuts required by these two formulations to solve the corresponding problems. As before, we have organized the information in a 4 ×4 grid of boxplox graphs. The reader can easily observe that sSEC always requires less number of cuts (blue boxes corresponding to SEC are always above orange ones corresponding to sSEC) showing that this formulation is more accurate than SEC: it reports smaller gaps (see Fig. 14 ) and needs less number of cuts. 8. Concluding remarks This paper has analyzed a novel version of the crossing postman problem with neighborhoods. We have shown that the problem can be cast within the framework of the family of mixed integer second order cone programming and several exact formulations are presented and computationally tested on an extensive testbed of instances. Additionally, we have presented a heuristic algorithm providing good quality solutions. It can be considered for large scale problems and also as a procedure to obtain initial solutions to initialize exact solvers handling our formulations. Computational results show that the problem is very hard and already for problems with 20 neighborhoods exact approaches fail to find an optimal solution within two hours of cpu time. This research opens up several research lines and extensions of the basic problem that can be included in the model. Among them 17 J. Puerto and C. Valverde European Journal of Operational Research xxx (xxxx) xxx ARTICLE IN PRESS JID: EOR [m5G; July 26, 2021;4:30 ] Fig. 14. Number of SEC added in the execution time. we mention finding better formulations or decomposition schemes that help in solving exactly larger instance sizes; and alternative heuristic algorithms that allow tackling large scale problems. Other extensions of the proposed models considered in this paper are the consideration of barriers that represent some buildings that the drone tour can not cross, conditions that control the displacement on the border of nonlinear neighborhoods like circles or problems that take into account the limited autonomy of drones requiring that the drone comes back to a depot to be recharged before to complete the route. Some of these topics will be the subject of a follow up paper. Acknowledgments This research has been partially supported by Spanish Ministry of Education and Science/FEDER grant number MTM2016-74983C02-(01–02), and projects FEDER-US-1256951, Junta de Andalucía P18-FR-1422, CEI-3-FQM331 and NetmeetData: Ayudas Fundación BBVA a equipos de investigación científica 2019. References Amorosi, L., Chiaraviglio, L., D’Andreagiovanni, F., & Blefari-Melazzi, N. (2018). Energy-efficient mission planning of UAVs for 5G coverage in rural zones. In Proceedings of the IEEE international conference on environmental engineering (ee) (pp. 1–9). https://doi.org/10.1109/EE1.2018.8385250 . Arkin, E. M., & Hassin, R. (1994). Approximation algorithms for the geometric covering salesman problem. Discrete Applied Mathematics, 55 (3), 197–218. https: //doi.org/10.1016/0166-218x(94)90 0 08-6 . Benders, J. F. (1962). Partitioning procedures for solving mixed-variables programming problems. Numerische Mathematik, 4 , 238–252. https://doi.org/10.1007/ BF01386316 . Blanco, V. (2019). Ordered p-median problems with neighbourhoods. Computational Optimization and Applications, 73 , 603–645. https://doi.org/10.1007/ s105890190 0 077-x . Blanco, V., Fernández, E., & Puerto, J. (2017). Minimum spanning trees with neighborhoods: Mathematical programming formulations and solution methods. European Journal of Operational Research, 262 (3), 863–878. https://doi.org/10.1016/ j.ejor.2017.04.023 . Blanco, V., & Puerto, J. (2021). On hub location problems in geographically flexible networks. International Transactions in Operational Research, n/a (n/a). https://doi. org/10.1111/itor.12993 . (2015). In Á. Corberán, & G. Laporte (Eds.), Arc routing: Problems, methods, and applications . Philadelphia (USA): Society for Industrial and Applied Mathematics. https://doi.org/10.1137/1.9781611973679 . 18 J. Puerto and C. Valverde European Journal of Operational Research xxx (xxxx) xxx ARTICLE IN PRESS JID: EOR [m5G; July 26, 2021;4:30 ] Díaz-Báñez, J. M., Mesa, J. A ., & Schöbel, A . (2004). Continuous location of dimensional structures. European Journal of Operational Research, 152 (1), 22–44. https://doi.org/10.1016/S0377-2217(02)00647-1 . Disser, Y., Mihalák, M., Montanar, S., & Widmayer, P. (2014). Rectilinear shortest path and rectilinear minimum spanning tree with neighborhoods: 8596 LNCS (pp. 208– 220). Springer Verlag. https://doi.org/10.1007/9783319091747 _ 18 . Edmonds, J. (2003). Combinatorial optimization eureka, you shrink! (pp. 11–26)). Springer . Garfinkel, R. , & Webb, J. (1999). On crossings, the crossing postman problem, and the rural postman problem. Networks: An International Journal, 34 (3), 173–180 . 10.1002/(SICI)1097-0037(199910)34:3%3C173::AID-NET1%3E3.0.CO;2-W Gentilini, I., Margot, F., & Shimada, K. (2013). The travelling salesman problem with neighbourhoods: Minlp solution. Optimization Methods and Software, 28 (2), 364– 378. https://doi.org/10.1080/10556788.2011.648932 . Geoffrion, A. M. (1972). Generalized Benders decomposition. Journal of Optimization Theory and Applications, 10 (4), 237–260. https://doi.org/10.10 07/BF0 0934810 . Gurobi Optimization, L. (2019). Gurobi optimizer reference manual, version 8.1.0. Knight, R. (2016). Drones deliver healthcare. http://insideunmannedsystems.com/ drones-deliver-healthcare/ . Lavars, N. (2015). Amazon to begin testing new delivery drones in the US. New Atlas, 13 . Lobo, M. S., Vandenberghe, L., Boyd, S., & Lebret, H. (1998). Applications of secondorder cone programming. Linear Algebra and its Applications, 284 (1), 193– 228. https://doi.org/10.1016/S0024-3795(98)10032-0 . International Linear Algebra Society (ILAS) Symposium on Fast Algorithms for Control, Signals and Image Processing Mallozzi, L., Puerto, J., & Rodríguez-Madrena, M. (2019). On location-allocation problems for dimensional facilities. Journal of Optimization Theory and Applications, 182 (2), 730–767. https://doi.org/10.1007/s1095701801470y . McCormick, G. P. (1976). Computability of global solutions to factorable nonconvex programs: Part I convex underestimating problems. Mathematical Programming, 10 , 147–175. https://doi.org/10.1007/BF01580665 . Miller, C. E., Tucker, A. W., & Zemlin, R. A. (1960). Integer programming formulation of traveling salesman problems. Journal of the ACM, 7 (4), 326329. https://doi. org/10.1145/321043.321046 . Mladenovi ´ c, N., & Hansen, P. (1997). Variable neighborhood search. Computers & Operations Research, 24 (11), 1097–1100. https://doi.org/10.1016/s0305-0548(97) 0 0 031-2 . Orloff, C. S. (1974). A fundamental problem in vehicle routing. Networks, 4 (1), 35–64. https://doi.org/10.10 02/net.3230 040105 . Otto, A., Agatz, N., Campbell, J., Golden, B., & Pesch, E. (2018). Optimization approaches for civil applications of unmanned aerial vehicles (UAVs) or aerial drones: A survey. Networks, 72 (4), 411–458. https://doi.org/10.1002/net.21818 . Pereira, V. (2018). Project: Metaheuristic-Local_Search-Variable_Neighborhood_ Search. https://github.com/Valdecy/Metaheuristic-Local _ Search-Variable _ Neighborhood _ Search . Puerto, J., & Valverde, C. (2021). Project: Instances for the crossing postman problem with neighborhoods (XPPN). https://github.com/z72vamac/examples _ xppn . Sawik, T. (2016). A note on the miller-Tucker-Zemlin model for the asymmetric traveling salesman problem. Bulletin of the Polish Academy of Sciences Technical Sciences, 64 (3), 517–520. https://doi.org/10.1515/bpasts20160057 . Schöbel, A. (2015). Location of dimensional facilities in a continuous space. In Location science (pp. 135–175). Springer International Publishing. https://doi.org/10. 1007/9783319131115 _ 7 . Velednitsky, M. (2017). Short combinatorial proof that the DFJ polytope is contained in the MTZ polytope for the asymmetric traveling salesman problem. Operations Research Letters, 45 (4), 323–324. https://doi.org/10.1016/j.orl.2017.04.010 . Yang, Y., Lin, M., Xu, J., & Xie, Y. (2007). Minimum spanning tree with neighborhoods: 4508 LNCS (pp. 306–316). Springer Verlag. https://doi.org/10.1007/ 9783540728702 _ 29 . Yuan, B., & Zhang, T. (2017). Towards solving TSPN with arbitrary neighborhoods: A hybrid solution. In Acalci (pp. 204–215). https://doi.org/10.1007/ 9783319516912 . 19