An Optimization Framework for Power Infrastructure Planning
AA N O PTIMIZATION F RAMEWORK FOR P OWER I NFRASTRUCTURE P LANNING
Nina Wiedemann
Institute for Operations ResearchETH ZurichSwitzerland [email protected]
David Adjiashvili
Institute for Operations ResearchETH ZurichSwitzerland [email protected]
February 2, 2021 A BSTRACT
The ubiquitous expansion and transformation of the energy supply system involves large-scale powerinfrastructure construction projects. In the view of investments of more than a million dollars perkilometre, planning authorities aim to minimise the resistances posed by multiple stakeholders.Mathematical optimisation research offers efficient algorithms to compute globally optimal routesbased on geographic input data.We propose a framework that utilizes a graph model where vertices represent possible locations oftransmission towers, and edges are placed according to the feasible distance between neighbouringtowers. In order to cope with the specific challenges arising in linear infrastructure layout, we firstintroduce a variant of the Bellman-Ford algorithm that efficiently computes the minimal-angle shortestpath. Secondly, an iterative procedure is proposed that yields a locally optimal path at considerablylower memory requirements and runtime. Third, we discuss and analyse methods to output k diversepath alternatives.Experiments on real data show that compared to previous work, our approach reduces the resistancesby more than 10% in feasible time, while at the same time offering much more flexibility andfunctionality. Our methods are demonstrated in a simple and intuitive graphical user interface, and anopen-source package (LION) is available at https://pypi.org/project/lion-sp . Power infrastructure layout is a challenge of significant impact on society. Besides ensuring consistent and robustenergy supply, it also plays a major role for the successful transition to renewable energies in view of climate change.The European Union, for example, identified necessary investments in the EU-transmission networks of around 50,000km overall length just from 2014 to 2024, of which 80% directly or indirectly target the integration of renewableenergies [13]. With an estimated cost of 1.5 million Euros per km [29], this leads to investments of approximately 75billion Euros over ten years. In light of this huge scale, it is more important than ever to optimize the planning process.But while planning authorities aim at minimal costs, environmentalists demand to consider the impact on conservationareas and bird mortality, and local residents oppose huge visible infrastructure. The planning of power lines has fittinglybeen called a ”field of tension” [5].Nevertheless, power infrastructure planning is oftentimes highly manual, mostly regarded as a communicative processbetween experts [5, 15]. In contrast, the optimal placement of power lines can be formulated as a simple combinatorialproblem, when provided with appropriate data extracted from geographic information systems (GIS). Possible computa-tional approaches for optimal layout of power lines have been discussed for a long time in the literature. In most casesthe area is rasterized and each cell is assigned a cost or resistance according to GIS data. For example, Popp et al. [22]presented very early work on optimal tower placement and even the selection of pylon heights based on the line sag and a r X i v : . [ m a t h . O C ] J a n PREPRINT - F
EBRUARY
2, 2021 l o w r es i s t a n ce G I S d a t a h i gh Figure 1:
An optimal transmission line planned with our framework. The resistances given with GIS-data are minimized togetherwith the angles along the route, yielding the globally optimal pylon placement. terrain. More recently, Bevanger et al. [5] proposed an environment-friendly power-line routing in Norway. A least-costpath (LCP) toolbox was published in the context of the same project [16]. Similarly, Monteiro et al. [20] and de Limaet al. [10] describe planning algorithms that consider the slope of the terrain and costs on direction change [20]. Othersdevelop decision support software that compute not only a single route, but multiple (pareto-optimal) alternatives, usingmulti-objective optimization techniques [3] or valley-finding algorithms applied on the cost surface [27].Despite the advances on optimizing power infrastructure layout, the approaches above suffer from a common drawback,namely that they only compute a shortest route through a regular grid, where each cell is connected to its 8-neighborhood.Transmission towers are only placed as a second step, yielding sub-optimal solutions: Crucially, it is disregarded thatthe cost of placing a power tower in a cell can differ significantly from the cost of traversing a cell with a cable. Weinstead suggest a graph model to minimize the resistance of the pylon placement, offering several advantages as detailedin Figure 2c.An additional challenge in power infrastructure planning is the avoidance of angles, since angles increase constructioncosts and electrical resistance. Work on angle-constrained shortest path computation in a graph model similar to ours isgiven by Santos et al. [26]. Building up on Piveteau et al. [21], they construct a line graph in a pre-computed corridorto model angle resistance. These works are important steps towards optimized pylon spotting, but leave room forimprovement with respect to computational efficiency. Contribution.
Here, we aim to fill this gap with an efficient framework for the globally optimal placement of pylons inpower transmission lines, offering the following functionality:• Graph modelling of transmission line planning with given (possibly different) cable and pylon resistances.• Provably efficient algorithms to compute the minimal-angle shortest path. The baseline approach, a linegraph [26, 21], has size O ( (cid:80) v ∈ V δ − v · δ + v ) for a graph with in-degree δ − v and out-degree δ + v at vertex v ∈ V .The Bellman Ford (BF) algorithm thus requires O ( | V | · (cid:80) v ∈ V δ − v · δ + v ) operations. In contrast, we adapt BFto operate on the normal graph with size O ( (cid:80) v ∈ V δ + v ) , and reduce the runtime to O (cid:0) | V | · (cid:80) v ( δ − v + δ + v ) · log( δ − v δ + v ) (cid:1) with a convex angle cost function.• An iterative procedure to achieve time and space feasibility with arbitrarily large project regions• Methods to compute a set of k paths that are short but also diverse among each other according to a suitablemetric. This is beneficial for planners as a decision support tool, to compare and analyse alternatives indifferent regions.We evaluate our methods on realistic power infrastructure projects in collaboration with Gilytics , a software companyworking on visualization and optimization software for linear infrastructure planning. We release an open-sourcepackage called LION (LInear Optimization Networks) providing all main algorithms. An exemplary output is shown inFigure 1. The term ”resistance” is used interchangeably with ”cost” here and thus does not refer to electrical resistance. Instead, it denotesthe summarized resistance of several geographic layers. PREPRINT - F
EBRUARY
2, 2021
The input to our framework is rasterized geographic data, here with a resolution of 10x10 meter cells. The ”resistanceraster” C , describing the costs per cell, is constructed from binary layers L , . . . , L r (Figure 1 left) and correspondingweights w , . . . , w r . Each layer indicates the occurrence of a geographic feature, e.g. L i [ x, y ] = 1 could indicatethat cell ( x, y ) is part of a forest. A feature L j that prohibits the placement of a transmission tower is integrated with w j = ∞ . We refer to previous works for approaches to decide how to weight and summarize these layers [5, 15]. Here,we assume that the resistance raster C is simply a weighted sum of the individual layers: C = r (cid:88) i =1 w i · L i (2.1) Similar to [26, 21], we model power infrastructure planning with a graph G = ( V, E ) with n := | V | verticesrepresenting transmission towers, and m := | E | directed or undirected edges, corresponding to the cable betweentowers. Figure 2a visualizes the graph layout: Each vertex is one cell in the raster, i.e. there exists a bijective function l : V → N × N , such that l ( v ) outputs the set of raster coordinates modelled by vertex v . Conversely, l − ( x, y ) is thevertex representing the cell at coordinates ( x, y ) . When not noted otherwise, s ∈ V will denote the source vertex and t ∈ V the target. Furthermore, there is an edge connecting two vertices u, v if and only if their respective pylons couldbe neighbors in a transmission line, which is the case if their distance is between two bounds d min and d max .With simplifying assumptions, the setting can be modeled in a directed acyclic graph (DAG) which significantlyaccelerates shortest path algorithms. Let (cid:126)a u,v := l ( u ) − l ( v ) be the vector from the cell of vertex u to cell of vertex v .We then place an edge between a pair of vertices if their distance is between d min and d max , but additionally requirethat (cid:126)a u,v diverges less than a certain maximum angle θ α from the straight line connection from source s to target t : E = (cid:110) ( u, v ) | d min ≤ (cid:107) (cid:126)a u,v (cid:107) ≤ d max ∧ cos (cid:0) (cid:126)a u,v · (cid:126)a s,t (cid:107) (cid:126)a u,v (cid:107)(cid:107) (cid:126)a s,t (cid:107) (cid:1) < θ α (cid:111) (2.2)This assumption is justified since in power infrastructure layout it is very unlikely that a tower is placed ”behind”another one in a transmission line from s to t . If θ α < π , then the graph is a DAG, because the neighbors of each vertexare the cells on the semicircular ring directed towards the target (Figure 2a). Note however that the algorithms presentedin the following sections are not restricted to DAGs. s t (cid:88) (cid:89) s t s ts t (a) Pylon spotting (ours) s t s t U V (b)
Line routing (two-stage approach often used in prior work)
Pylon spotting (ours) Line routing (prior work)
Global optimizationof pylon positions Two stage approachForbidden regions canbe traversed by cable Forbidden areascannot be crossedAllows to differentiateresistances forpylons and cable No differentiationpossibleWill find feasiblepylon spotting(if existent) Placing pylons alongthe precomputed routemight be impossibleDirect angleminimization Angles along route (cid:54) = angles at pylons (c) Advantages of our approach
Figure 2:
Comparison of our graph model (a) to previous work (b). In a) each raster cell is a vertex, and edges are placed betweeneach pair of vertices with distance between d min and d max . Only edges with sufficiently low angle to the straight line from s to t areallowed, thereby yielding a directed acyclic graph (DAG). In previous approaches (b), first the optimal route through a regular grid isfound, and only in the second step transmission towers are placed along the route. The drawbacks of this line-routing compared todirect pylon spotting are listed in c). In contrast to two-stage line routing approaches (Figure 2b) used in previous work [20, 30, 31, 18], a pylon-spotting approach allows a) to differentiate pylon and cable costs, b) to find a globally optimal pylon spotting and c) to minimize3
PREPRINT - F
EBRUARY
2, 2021the angles at the pylons. In Figure 2c both approaches are contrasted in more detail. Pylon and cable resistances canbe computed independently with different weights, such that C p [ x, y ] is the cost of placing a pylon in cell ( x, y ) and C c [ x, y ] denotes the cost to traverse ( x, y ) with a cable. In order to translate C p and C c into edge costs in a graph, thecable cost must be summed over several cells, since the cable resembles a straight line above a discrete raster of cells(black arrow in Figure 2a). To approximate this line, we utilize a well-known line drawing algorithm from computergraphics, the Bresenham line [6]. The computation of the Bresenham line M e for an edge e = ( u, v ) yields a setof coordinates corresponding to the raster cells on the line, formally M e = { l ( u ) , ( x , y ) , . . . , ( x d − , y d − ) , l ( v ) } .Pylon and cable costs are then combined with a weight w c that describes the importance of the cable cost (i.e. theaverage resistance of the cells in M e ) in contrast to the pylon cost: c ( e ) = C p [ l ( u )] + C p [ l ( v )]2 + w c · d (cid:88) ( x,y ) ∈ M e C c [ x, y ] (2.3) We use the Bellman-Ford (BF) [4, 14] algorithm to compute the shortest (minimal resistance) path, since it is morespace-efficient than Dijkstra and the runtime can be reduced if an upper bound on the number of vertices on the shortestpath is known. While the BF algorithm will minimize cable and pylon resistances according to Equation 2.3, it isnon-trivial to regard di-edge costs such as angle costs. The latter problem was formalized as a special case of the
Quadratic Shortest Path Problem (QSPP) [24], named the Adjacent QSPP. In AQSPP the objective is to minimize aninteraction cost of adjacent edges along with the normal edge costs. Even though AQSPP was shown to be NP-hard whenthe graph is not acyclic [25, 17], an auxiliary graph structure, namely the line graph or dual graph , has been exploitedin applications to deal with di-edge costs, including power line [26, 21] and other infrastructure planning [7, 32].However, a major drawback of the line graph is its size: A directed graph with δ − v incoming edges and δ + v outgoing edgesat vertex v results in (cid:80) v δ − v · δ + v edges in the line graph. Less memory-consuming representations have been proposedin work on turn constraints in road networks [11], where Dijkstra is modified to consider turn-tables associated witheach vertex. Here, we propose a variant of the Bellman-Ford algorithm named MABF, that computes the minimal-angle walk (or path under certain conditions) in the original graph. In addition to the edge weight given by the cost function c ( e ) : E −→ R in Equation 2.3, let α : [0 , −→ R + be anangle cost function. We aim to compute a distance map D and a predecessor map T from which the shortest path canbe reconstructed. Specifically, D [ e ] describes the accumulated costs for the (current) shortest path from the source s tothe edge e , including angle costs. D is initialized to infinity for all edges, except for the outgoing edges e = ( s, − ) ofthe source, which are instead set to their respective edge cost, formally D [ e ] = c ( e ) ∀ e = ( s, − ) . u v... (a) Update step in minimal-angle BF vust w (b)
Cycle on optimal walk
Figure 3:
Minimal-angle shortest path computation: As shown in a), the minimum over the incoming edges f , . . . , f d must betaken to find the optimal predecessor of edge e with respect to the distance D and angle cost α . For non-concave α , the algorithmcan yield cyclic paths as in b), since high angle costs (red) can be avoided by using more but smaller angles along the cycle (greenand orange). Algorithm 1 outlines the proposed procedure to find the angle-penalized least cost walk from s to t . Note that in thegeneral case the route can contain cycles because thereby sharp turns (high angle costs) are avoided (see Figure 3b). Wefirst describe the algorithm and then show that the result will be a path when α is concave. We assume an upper bound p on the number of edges in the shortest path. Note that p < n , but due to the raster-layout of the power infrastructuregraph, it is usually safe to set p to a value much smaller than n . BF computes the shortest path tree by updating the4 PREPRINT - F
EBRUARY
2, 2021distances to all edges for p iterations. Figure 3a visualizes the new update step applied to edge e = ( u, v ) : Its distancefrom s is given by the minimum over the distances of its possibly preceding edges, together with the respective anglecost: D [ e ] = min f ∈ E, f =( − ,u ) D [ f ] + c [ e ] + α ( (cid:94) ( f, e )) (3.1) Algorithm 1:
Minimal-angle Bellman-Ford algorithmInput: G , s , c , αD [ e ] = ∞ ∀ e ∈ ED [( s, − )] = c (( s, − )) T [ e ] = e ∀ e ∈ E for p iterations dofor i ← to m do ( u, v ) ← e i D [ e i ] = min f ∈ E, f =( − ,u ) D [ f ] + c [ e i ] + α ( (cid:94) ( f, e i )) T [ e i ] = argmin f ∈ E, f =( − ,u ) D [ f ] + c [ e i ] + α ( (cid:94) ( f, e i )) endend In the following, it is shown that this variant of the Bellman-Ford algorithm computes the correct distances.
Theorem 1.
Let e t = ( − , t ) be any edge in G . Furthermore, let P ∗ be a least cost walk of length | P ∗ | ≤ p from s to t ,and denote the edges on P ∗ as P ∗ = { e , ..., e t } .After p update iterations in Algorithm 1, it holds that D [ e t ] = c ( e ) + t (cid:88) i =2 c ( e i ) + α ( (cid:94) ( e i − , e i )) := c ( P ∗ ) and c ( P ∗ ) is the cost of the shortest minimal-angle walk.Proof. Proof by induction on p . The base case p = 1 is given by the initialization D [( s, − )] = c (( s, − )) . As theinduction step, consider the least cost walk P ∗ from s to e t = ( u, t ) such that ( c ( P ∗ ) is minimal. Let f ∗ be the edgepreceding e t on P ∗ , so f ∗ is of the form ( − , u ) . Then P ∗ \ e t is a shortest walk from s to f ∗ of length p − . Byinduction, after p − iterations we have D [ f ∗ ] = c ( P ∗ \ e t ) . More general, for all incoming edges of u , defined as f j : f j ∈ E, f j = ( − , u ) , D [ f j ] contains the total cost of the shortest walk P j from s to f j satisfying | P j | ≤ p − .Then the distance of e t is updated according to Equation 3.1: D [ e t ] = min j D [ f j ] + c [ e t ] + α ( (cid:94) ( f j , e t ))= D [ f ∗ ] + c [ e t ] + α ( (cid:94) ( f ∗ , e t ))= c ( P ∗ \ e t ) + c [ e t ] + α ( (cid:94) ( f ∗ , e t ))= c ( P ∗ ) Thus, f ∗ is correctly assigned as the predecessor of e t . As shown in Figure 3b, the route with least angle resistances might contain cycles, since strong angles (red) are avoidedwith multiple small angles (green, orange). The following lemma proves a sufficient condition that guarantees that theoutput is a path . Lemma 1.
If the angle cost function α : [0 , −→ R + is concave, algorithm 1 returns a shortest path.Proof. According to the Closed-Path theorem formulated by Abelson and DiSessa [1], ”a total turning along any closedpath is an integer multiple of ◦ ”, which applies to a cycle in a digraph. However, the angle at the entry vertex5 PREPRINT - F
EBRUARY
2, 2021is not included in our case, leaving at least ◦ for the cycle. Let σ , . . . , σ n denote all angles encountered on thecycle. Then ∀ i : σ i ∈ [0 , and (cid:80) ni =1 σ i ≥ due to the Closed-Path theorem. It is further well-known that for aconcave function α it holds that (cid:80) i α ( x i ) ≥ α ( (cid:80) i x i ) , and therefore it follows n (cid:88) i =1 α ( σ i ) ≥ α (cid:16) n (cid:88) i =1 σ i (cid:17) ≥ α (180) . We conclude that the angle costs accumulating on the cycle will be at least as large as the angle cost that would be partof the direct path, which is at most α (180 ◦ ) . The above results show the correctness of algorithm 1 to compute the minimal-angle shortest walk for general andthe shortest path for concave angle cost functions. While it improves on the space requirements in contrast to a linegraph construction, namely O ( m ) in contrast to O ( (cid:80) v δ − v · δ + v ) , it has the same runtime: In each update of an edge, theminimum over the incoming edges is taken, leading to a runtime of O (cid:0) p · (cid:80) v δ − v · δ + v (cid:1) , corresponding to the runtime ofnormal BF applied on the line graph. We however claim that the runtime can be improved significantly, if the angle costfunction is a discrete (step) function or convex. Lemma 2. If α : [0 , −→ S is a step function with | S | possible angle costs ( | S | < δ − v , δ + v ), the shortest walk canbe computed in time O (cid:0) p · (cid:80) v ( | S | · δ − v + δ + v ) · log( | S | δ − v δ + v ) (cid:1) . Lemma 3.
For convex, monotonically increasing α the shortest walk can be computed in time O (cid:0) p · (cid:80) v ( δ − v + δ + v ) · log( δ − v δ + v ) (cid:1) . The proofs for Lemma 2 and Lemma 3 are found in Appendix 2 and B respectively. Both algorithms are implementedin the LION package and are used automatically based on the size of the instance and the selected angle cost function.To summarize, we provide space- and time efficient shortest path algorithms that minimize angles along the path andomit the need for a line graph construction.
Despite space-efficient algorithms, the graphs can still become extremely large, in the order of billions of edges. Inparticular, m increases quartic with the resolution: Not only the number of vertices grows quadratically, but also thenumber of neighbors per vertex. We take advantage of this observation with a multi-scale approach that narrows downthe region of interest iteratively.Initially, the path is computed on a coarse downsampled version of the instance. Then the resolution is increased overseveral steps while simultaneously restricting the project region to a corridor around the previously computed path.Let ( r , . . . , r q ) be a list of downsampling factors, and let ϑ be an upper bound on the number of edges (i.e. a memorylimit). Figure 4 provides a visual explanation: The resolution of the instance is decreased by a factor of r and theoptimal path P is found at this scale. For the next iteration, the project region is restricted to a corridor around P Figure 4:
Iterative multi-scale approach: The instance is downsampled by r and the optimal path P is computed (left). For the nextiteration, the region is restricted to a corridor of width w around P (dotted line left), and again the best path in the corridor withresolution r is yielded. The output is the best path through the final region (resolution r in corridor with size w ). PREPRINT - F
EBRUARY
2, 2021(white dotted line). In this corridor, the optimal path is computed with resolution corresponding to r . These steps arerepeated until reaching the highest resolution (usually r q = 1 ). In each step i , the width w i of the corridor is set to themaximum feasible width such that the number of edges remains below ϑ in the next iteration (where the instance isdownsampled by r i +1 ). Importantly, this multi-scale approach can also be combined with the computation of k diverseshortest paths explained in the following section: The corridor for the next iteration is in this case a buffer around all k paths. In general, the multi-scale procedure will yield only a locally optimal path, but with a sufficiently homogeneousinstance, good results can be achieved in significantly lower time and space. Since the single optimal path is solely based on the geographic input data (hard criteria), it might not comply with softcriteria coming up in the analysis of experienced planners. To trade off such criteria in discussions between plannersand involved stakeholders, it is of interest to propose a set of multiple short paths for the sake of comparison.A very efficient method to compute k shortest paths is given with Eppstein’s algorithm [12] if paths are allowed to beloopy. Here, it is first shown that Eppstein’s algorithm can be easily adapted to return the k minimal-angle shortestpaths. As defined above, let p be an upper bound on the path length. Lemma 4.
The k shortest paths can be computed in O ( p · (cid:80) v δ − v · δ + v + m log m + kp ) Proof.
The proof is given by the following adaptation of Eppstein’s algorithm (EA): As in EA, not only the distance D s and predecessor map T s for the shortest path tree rooted in s , but also the tree rooted in t is computed with reversededge directions, yielding D t , T t . With algorithm 1 this takes O (cid:0) p · (cid:80) v δ − v · δ + v (cid:1) operations. Following EA, D s and D t are summed up point-wise and the edge costs c ( e ) are subtracted since they appear twice. The resulting map S , with S [ e ] = D s [ e ] + D t [ e ] − c ( e ) , now describes the distance of the shortest path through e . Then S is sorted O ( m log( m )) ,and the paths are reconstructed for the edges with the k lowest values in S , which can be done in O ( kp ) using T s and T t .Enumerating the k shortest paths often leads to highly similar paths. Thus, the problem of diverse path finding hasbeen introduced by Liu et al. [19] as the ”Top- k Shortest Paths with Diversity (KSPD)”, which has been shown tobe NP-hard [9]. Both Liu et al. [19] and Chondrogiannis et al. [8] provide greedy methods to find k paths that aresufficiently diverse with respect to a threshold θ of a suitable path diversity metric. Akgün et al. [2] on the other handconsider the reversed formulation, namely to find the k most diverse paths with a threshold on the path costs, therebyturning the problem into a p-dispersion problem. Here, we were first interested to understand for which metrics it iseven possible to find one diverse path efficiently. Definition 1.
Given a graph G = ( V, E ) , a path P and a metric M , we define the diverse path selection problem asthe problem to find an alternative path ˆ P where the distance exceeds a similarity threshold θ , i.e. d M ( P, ˆ P ) > θ . We first consider the Yau-Hausdorff distance d Y , which is the maximum Euclidean distance between two sets P , Q ,and was shown to be a metric by Tian et al. [28]. d Y ( P, Q ) = max { max p ∈ P min q ∈ Q (cid:107) p − q (cid:107) , max q ∈ Q min p ∈ P (cid:107) p − q (cid:107)} (5.1) Lemma 5.
Diverse path selection according to the Yau-Hausdorff distance is decidable in polynomial time.Proof.
After computing both shortest path trees, one pass over all edges is required to determine which vertices havesufficient Euclidean distance to the previous path(s) P prev . Note that in order to solve the diverse path selection problem, any diverse path is sufficient, but with Eppstein’s algorithm we can also find the shortest diverse path by finding thesufficiently diverse path with minimal S [ e ] : e ∗ = argmin e =( u,v ): min w ∈ Pprev (cid:107) u − w (cid:107) >θ S [ e ] . If e ∗ = ( u ∗ , v ∗ ) exists and has non-infinity path distance, the shortest path P ∗ through e ∗ can be computed from T s and T t . The Yau-Hausdorff distance of P ∗ to P prev is the maximum distance of all vertices of P ∗ , and thus at least aslarge as the distance of u ∗ to P : d Y ( P ∗ , P prev ) ≥ max p ∈ P ∗ min q ∈ P prev (cid:107) p − q (cid:107) ≥ min q ∈ P prev (cid:107) u ∗ − q (cid:107) Conversely, if no such edge exists then clearly there is no path with sufficient Yau-Hausdorff distance, because themaximum distance must be attained at some vertex. 7
PREPRINT - F
EBRUARY
2, 2021
Figure 5:
NP-hardness of finding the shortest path with a lower bound on the mean Euclidean distance or the Jaccard distance. If adiverse path exists, G contains a Hamilton path.
Lemma 5 thus suggests an efficient greedy algorithm to compute k diverse shortest paths with respect to d Y . Theruntime is given by computing the shortest path trees and subsequently computing the path distance ( O ( kp ) ) for all n vertices for k times, resulting in O ( mpd + k np ) .On the other hand, despite the proposed methods and promising computational results in [9, 19], we find that diversepath selection is NP-hard for the Jaccard distance metric d J and the mean Euclidean distance d M , d M ( P, Q ) = max { | P | (cid:88) p ∈ P min q ∈ Q (cid:107) p − q (cid:107) , | Q | (cid:88) q ∈ Q min p ∈ P (cid:107) p − q (cid:107)} (5.2) Lemma 6.
With the Jaccard distance or the average Euclidean distance metric, the diverse path selection problem isNP-hard.Proof.
Given a graph G = ( V, E ) and n = | V | , we construct an auxiliary graph H = ( V (cid:48) , E (cid:48) ) as in Figure 5: V (cid:48) = V ∪ { v (cid:48) i | i ∈ [1 ..n ] } E (cid:48) = E ∪ (cid:110) ( v (cid:48) i , v (cid:48) i +1 ) | i ∈ [1 ..n − (cid:111) ∪ (cid:110) ( v i , v (cid:48) (cid:4) n (cid:5) ) | ∀ i ∈ [1 ..n ] (cid:111) ∪ (cid:110) ( v i , v (cid:48) (cid:6) n +12 (cid:7) ) | ∀ i ∈ [1 ..n ] (cid:111) G can be any graph, but here we arrange all its vertices in a straight line in an Euclidean space, and place v (cid:48) . . . v (cid:48) n in aline parallel to v . . . v n such that d M ( v i , v (cid:48) i ) = 2 ∀ i .We claim that the Hamiltonian problem can be decided for G if it can be decided whether there exists a diversealternative to P = ( v (cid:48) , v (cid:48) , . . . , v (cid:48) n ) (horizontal green line in Figure 5). The threshold for the minimum Jaccard distanceis set to θ = 0 . . Observe that all vertices of P will also be part of the alternative path, and only the edges from v (cid:48) (cid:6) n +12 (cid:7) and v (cid:48) (cid:4) n (cid:5) allow the path to traverse other vertices than P . The Jaccard distance corresponds to the inverted IoU, so θ = 0 . requires the path to use at least | P | new vertices. But since | P | = n , the only possible alternative is a path thatvisits all vertices of G , i.e. a Hamiltonian path in G . If no alternative is found, it is implied that there is no Hamiltonpath because the connections from v (cid:48) (cid:6) n +12 (cid:7) and v (cid:48) (cid:4) n (cid:5) to all vertices of G ensure that it would be found if existent.Similarly, with respect to the mean-Eucledian distance d M we set θ = 1 . If k vertices of G are used ( k ≤ n ), then d ( P, P ∗ ) = n · k · n + k = ⇒ ( d ( P, P ∗ ) = 1 ⇐⇒ k = n ) Thus, a diverse alternative exists iff a Hamiltonian path exists.Nevertheless, one can still utilize Eppstein’s algorithm to enumerate paths sorted by distance, and then select the onesgreedily that have sufficient average Euclidean or Jaccard distance to the previously chosen paths, similarly to theFindKSP algorithm by Liu et al. [19]. We therefore implemented and compared five computational methods:• find-ksp-max : greedy path selection according to d Y , guaranteed to give the correct result in polynomial time8 PREPRINT - F
EBRUARY
2, 2021• find-ksp-mean : greedy path selection according to d M • greedy-set : greedy path selection according to d J • k-dispersion : 2-approximation [23] for p-dispersion problem, i.e. find most diverse set of paths given athreshold on the cost (similar to [2])• corridor-penalizing : Soft variant of find-ksp-max , where the vertices that are close to the previous paths areonly penalized instead of excludedExamples of five diverse minimal-angle and least resistance paths, computed with find-ksp-max , are shown in Figure 6b. In collaboration with the software company Gilytics we evaluate our methods on two transmission line planningprojects, labeled instance 1 and instance 2 in the following.
Instance 1 is an exemplary project in Austria with GISdata publicly available (resistance raster shown in Figure 6b). Instance 2 is a former project in Belgium for whichdata was extracted from OpenStreetMap . The resistances and an optimal path is shown in Figure 1. For both instances,realistic resistance values were decided together with Gilytics. Further details on the instances are given in Table 1.
Instance Rasterwidth Rasterheight Project region(not forbidden) d_min d_max neighboursper node1
25 35 967
15 25 632
Table 1:
Instance properties. All numbers refer to a raster of 10x10m cells.
All experiments were executed on a compute node equipped with two 64-core AMD EPYC 7742 processors (2.25 GHznominal, 3.4 GHz peak) and 512 GB of DDR4 memory, clocked at 3200 MHz. In consultation with the collaborator itwas decided to assume a directed acyclic graph (DAG) according to the criterion in Equation 2.3. Note that all proposedalgorithms are also applicable for a general graph.
While our general modelling approach is similar to [26, 21], our framework allows to process much larger instancesdue to efficient algorithms omitting the need for larger graph representations such as a line graph. Specifically, fora single path computation, the required memory mainly stems from the storage of distance and predecessor maps,together comprising m elements. For finding k paths on the other hand, m elements are kept in RAM since twoshortest path trees are computed. With respect to actual memory, distances and predecessors are represented as integeror floating points of bytes, resulting in m · m bytes or m bytes respectively. In practice, m depends on theproject region and d min and d max . Here, processing instance 2 with m ≈ · edges (see Table 1) requiresaround . GB, but in our framework only edges reachable from the source are stored, such that the optimal path andeven the k shortest paths can be computed with a standard 8GB RAM. In contrast, if instance 2 was processed with aline graph [26, 21], the number of edges and thus also the required memory would increase by a factor of 632. Thisdemonstrates the large gain in efficiency achieved with our framework. On top of that, any conceivable instance canbe processed with our multi-scale approach. Figure 6b depicts k shortest paths on instance 1 , computed with themulti-scale approach, which would otherwise require · · ≈ GB memory
For instance 1 , we extracted the pylon positions of an existing transmission line, denoted by Q in the following. InFigure 6b, Q is shown in contrast to five paths P , . . . , P computed with our framework. Although we naturally do notknow which geographic features were considered in the planning of Q , we can analyze Q with respect to our data, as it PREPRINT - F
EBRUARY
2, 2021 (a)
Comparison of resistance of the optimal path computed with LION and the ex-isting power line. Clearly, the existing power line is traversing areas of higher resis-tances than our path. optimal pathalternative 1alternative 2alternative 3alternative 4existing power line (b)
K diverse shortest paths for instance 1 (cropped) are shown, computed withthe find-ksp-max method. Dark areas indicate low resistances, brighter areashigher resistance, and white areas are forbidden. High resistance areas areavoided by our paths in contrast to the existing power line
Figure 6:
Qualitative and quantitative comparison of our results with existing infrastructure is safe to assume that similar resistances were regarded when planning Q . Figure 6a demonstrates that Q is crossingmany high resistance areas such as bird sanctuaries or even areas that are forbidden according to our data (buildings,nature reserves), while P tends to follow roads and suitable regions around forest. Even when omitting costs forpylons that traverse forbidden regions, we find that the geographic resistances are reduced by 25%, together with an 8%reduction of angles, rendering P a strictly better route. As explained above (Figure 2), most related works [20, 30, 31, 18] compute a route through the resistance raster andplace transmission towers as a second step (called line-routing in the following), instead of computing a globally optimal pylon spotting . In Table 2, we compare the outputs of our pylon spotting to the line-routing approach. Specifically,we simply apply our algorithm on the graph with d min = 1 , d max = 1 . , i.e. a grid-structured graph where each cellis connected to its 8-neighborhood. When the angle weight was greater than zero, we reduced also the angles of theline-route with our algorithm to provide a fair comparison of the angle costs. Based on this route, the pylons are spotted,where we even minimize the resistances again with our algorithm, but limiting the project region to locations along theroute (Figure 2b). Instance (resolution)
Angleweight Approach Time(seconds) m (|E|)(mio) Anglecost Resistance Distance(in m)Line-routing 34.1
Line-routing 35.6
Line-routing 33.1
Inst. 1 (20 m)
Line-routing 23.8
Line-routing 21.9
Line-routing 20.0
Inst. 2 (10 m)
Table 2:
Comparison of pylon-spotting to line routing (best value marked bold). ”Distance” is the average distance between thepylons of our approach compared to the baseline. Although our approach leads to much larger graphs, the optimal path can still becomputed in less than 10 minutes and in all cases leads to significantly lower angle and resistance costs. PREPRINT - F
EBRUARY
2, 2021Although our graphs are larger than grid-graphs by a few orders of magnitude, the instances can still be processed inless than 10 minutes . The increase in runtime is, however, rewarded by the achieved reduction in costs: Resistancesdecrease by 10-20% while the angles are reduced by more than 50%. A major drawback of line-routing is also thatminimizing the angle along the route does not necessarily reduce the angles at the pylons, but can in fact result in largerangles despite penalizing angles more (see instance 2 angle weight 0.1). Furthermore, the routes are substantiallydifferent (see ”distance” in Table 2), such that a local search around the route does not suffice to determine the optimalpath. We explored in different configurations how much the space and time could be reduced by our itera-tive multi-scale method (section 4), and at what cost. In Figure 7 the performances for instance 2 are reported for the scale lists [2 , , [3 , and [3 , , , meaning for example that the instance is down-sampled to 30x30m resolution at first, then to 20x20m and last to 10m. The bars show the relative increase or decrease of runtime, space and resistance, always compared to the direct computation ( [1] ). [1] [2, 1] [3, 1] [3, 2, 1] List of scales in multi-scale processing P e r c e n t o f o p t i m a l p a t h Runtime (s) Cost Max. number of edges
Figure 7:
Evaluation of iterative multi-scale methods to reduce time and spacerequirements. The outputs were computed for three instances with four differentlist of scales ([1] refers to a direct optimization in one iteration). Space andruntime decrease significantly for all instances when using the iterative procedure,and the resistances increase only marginally as desired.
Importantly, the path resistances only in-crease by less than 5%, while in some cases(e.g. [2 , ), the optimal path is found. Im-portantly, space time complexity is reduceddramatically due to the quartic dependenceof the number of edges on the raster resolu-tion. Specifically, down-scaling by factors and decreases the number of edges byfactors and , respectively. The overallruntime in all cases was below 13% of theoriginal runtime. Furthermore, we observedthat having three scales, instead of two (e.g. [3 , , vs. [3 , ) increases runtime, but alsoleads to slightly lower path resistances, indi-cating that more steps improve the ability toapproximate the optimal path. These results,summarized in Figure 7, demonstrate the po-tential of the multi-scale approach to findclose-to-optimal paths at low computationalcosts. In section 5 several methods to compute the k diverse shortest paths were discussed. Here we apply these methods on instance 2 and vary the diversity threshold, observing the ability of each method to achieves high diversity at low pathcosts. Each point in Figure 8 corresponds to a set of 5 paths, and the diversity of these paths is plotted against the sumof their respective path costs. Optimally the points should be located in the lower right corner, such that they minimizecosts but maximize diversity in the set of paths.As expected, find-ksp-max and corridor-penalizing obtain best results with respect to the Yau-Hausdorff metric,while greedy-set and k-dispersion algorithms achieve high Jaccard distance at low costs. In short, guarantees are onlygiven with find-ksp-max , but other algorithms such as greedy-set are shown to be effective and should be employeddependent on the diversity metric. instance 1 was down sampled to 20m resolution to allow for optimal path computation without the multi-scale approach PREPRINT - F
EBRUARY
2, 2021
Jaccard distance s u m o f c o s t s Yau-Hausdorff distance (in raster cells)
Algorithm (metric) corridor-penalizing (Y)find-ksp-max (Y)find-ksp-mean (M)greedy-set (J)k-dispersion (Y)k-dispersion (J)
Mean Euclidean distance (in raster cells)
Figure 8:
Diversity-cost trade-off of various algorithms ( instance 2 ). Each point refers to a set of five paths that were computed witha certain algorithm and metric (Y - Yau-Hausdorff, J - Jaccard, M - Mean Euclidean distance). The greedy-set method minimizesthe Jaccard distance, while corridor-penalizing or find-ksp-max perform best for Euclidean metrics. The layout of power transmission lines is a contentious topic with many stakeholders involved, calling for softwareand mathematical methods to support the planning process. In contrast to previous work, we presented an efficientframework for the global optimization of transmission line layout, for computing (angle-) optimal pylon spotting. Onthe one hand, we view our methods as a mathematically sound way to propose a route that is optimal with respect to theinput parameters, while on the other hand it can also be seen as a decision support system for planners, with its abilityto output diverse shortest paths for comparison.While our experiments show the system’s ability to reduce resistances and angle costs compared to previous work,further work should aim at accelerating the proposed algorithms with approximation algorithms or parallelization. Manyfurther challenges remain, including for example the optimization of pylon heights in more complex (non-flat) terrain.
Figure 9:
User interface for power infrastructure planning
Alongside a
PyPi package lion-sp , we also re-lease a small graphical user interface , shown inFigure 9, demonstrating how a decision supporttool can be created by adding appropriate UI/UXcomponents. Moreover, the proposed method-ology be be extended beyond transmission lineplanning. Other applications such as the lay-out of pipes, underground transmission lines,or highways could benefit from the same opti-mization framework. We thus hope to pave theroad towards software-assisted linear infrastruc-ture planning, backed by proven mathematicalmethods, for easier collaboration between localresidents, planning agencies, environmentalistsand construction authorities. Acknowledgments
We are very grateful for the fruitful collaboration with Gilytics AG, and want to deeply thank them for their support incollecting the relevant data and understanding the background of the problem, and for their extremely helpful feedbackon our methods. Special thanks to Stefano Grassi for his contribution.Finally, we also want to thank Rico Zenklusen for his support, as well as Jordi Castells for his great help on improvingthe LION package. Installation instructions available on GitHub https://github.com/NinaWie/PowerPlanner PREPRINT - F
EBRUARY
2, 2021
References [1] H. Abelson and A. A. DiSessa.
Turtle geometry: The computer as a medium for exploring mathematics . MITpress, 1986.[2] V. Akgün, E. Erkut, and R. Batta. On finding dissimilar paths.
European Journal of Operational Research , 121(2):232–246, 2000.[3] D. Bachmann, F. Bökler, J. Kopec, K. Popp, B. Schwarze, and F. Weichert. Multi-objective optimisation basedplanning of power-line grid expansions.
ISPRS International Journal of Geo-Information , 7(7):258, 2018.[4] R. Bellman. On a routing problem.
Quarterly of applied mathematics , 16(1):87–90, 1958.[5] K. M. Bevanger, G. Bartzke, H. Brøseth, et al. Optimal design and routing of power lines; ecological, technicaland economic perspectives (optipol). progress report 2010.
NINA rapport , 2010.[6] J. E. Bresenham. Algorithm for computer control of a digital plotter.
IBM Systems journal , 4(1):25–30, 1965.[7] T. Caldwell. On finding minimum routes in a network with turn penalties.
Communications of the ACM , 4(2):107–108, 1961.[8] T. Chondrogiannis, P. Bouros, J. Gamper, and U. Leser. Alternative routing: k-shortest paths with limited overlap.In
Proceedings of the 23rd SIGSPATIAL International Conference on Advances in Geographic InformationSystems , pages 1–4, 2015.[9] T. Chondrogiannis, P. Bouros, J. Gamper, U. Leser, and D. B. Blumenthal. Finding k-dissimilar paths withminimum collective length. In
Proceedings of the 26th ACM SIGSPATIAL International Conference on Advancesin Geographic Information Systems , pages 404–407, 2018.[10] R. M. de Lima, R. Osis, A. R. de Queiroz, and A. H. M. Santos. Least-cost path analysis and multi-criteriaassessment for routing electricity transmission lines.
IET Generation, Transmission & Distribution , 10(16):4222–4230, 2016.[11] D. Delling, A. V. Goldberg, T. Pajor, and R. F. Werneck. Customizable route planning. In
International Symposiumon Experimental Algorithms , pages 376–387. Springer, 2011.[12] D. Eppstein. Finding the k shortest paths.
SIAM Journal on computing , 28(2):652–673, 1998.[13] O. Feix, R. Obermann, M. Strecker, and R. König. Netzentwicklungsplan strom 2014: Erster entwurf derübertragungsnetzbetreiber. , 2014.[14] L. R. Ford Jr. Network flow theory. Technical report, Rand Corp Santa Monica Ca, 1956.[15] T. Grossardt, K. Bailey, and J. Brumm. Analytic minimum impedance surface: Geographic information system-based corridor planning methodology.
Transportation Research Record , 1768(1):224–232, 2001.[16] F. Hanssen, R. May, J. Van Dijk, et al. Spatial multi-criteria decision analysis (smcda) toolbox for consensus-basedsiting of powerlines and wind-power plants (consite).
NINA Rapport , 2018.[17] H. Hu and R. Sotirov. Special cases of the quadratic shortest path problem.
Journal of Combinatorial Optimization ,35(3):754–777, 2018.[18] Y. Li, Q. Yang, W. Sima, J. Li, and T. Yuan. Optimization of transmission-line route based on lightning incidencereported by the lightning location system.
IEEE transactions on power delivery , 28(3):1460–1468, 2013.[19] H. Liu, C. Jin, B. Yang, and A. Zhou. Finding top-k shortest paths with diversity.
IEEE Transactions on Knowledgeand Data Engineering , 30(3):488–502, 2017.[20] C. Monteiro, I. J. Ramírez-Rosado, V. Miranda, P. J. Zorzano-Santamaría, E. García-Garrido, and L. A. Fernández-Jiménez. Gis spatial analysis applied to electric line routing optimization.
IEEE transactions on Power Delivery ,20(2):934–942, 2005.[21] N. Piveteau, J. Schito, M. Raubal, and R. Weibel. A novel approach to the routing problem of overhead transmissionlines. Master’s thesis, Geographisches Institut der Universität Zürich, 2017.[22] R. Popp, C. Dabekis, and F. Fullerton. Electronic computer program permits optimized spotting of electrictransmission towers.
IEEE Transactions on Power Apparatus and Systems , 82(66):360–365, 1963.[23] S. Ravi, D. J. Rosenkrantz, and G. K. Tayi. Facility dispersion problems: Heuristics and special cases. In
Workshopon Algorithms and Data Structures , pages 355–366. Springer, 1991.[24] B. Rostami, F. Malucelli, D. Frey, and C. Buchheim. On the quadratic shortest path problem. In
InternationalSymposium on Experimental Algorithms , pages 379–390. Springer, 2015.13
PREPRINT - F
EBRUARY
2, 2021[25] B. Rostami, A. Chassein, M. Hopf, D. Frey, C. Buchheim, F. Malucelli, and M. Goerigk. The quadratic shortestpath problem: complexity, approximability, and solution methods.
European Journal of Operational Research ,268(2):473–485, 2018.[26] A. H. M. Santos, R. M. de Lima, C. R. S. Pereira, et al. Optimizing routing and tower spotting of electricitytransmission lines: An integration of geographical data and engineering aspects into decision-making.
ElectricPower Systems Research , 176:105953, 2019.[27] J. Schito, D. Moncecchi, and M. Raubal. Determining transmission line path alternatives using a valley-findingalgorithm.
Computers, Environment and Urban Systems , 86:101571, 2021.[28] K. Tian, X. Yang, Q. Kong, C. Yin, R. L. He, and S. S.-T. Yau. Two dimensional yau-hausdorff distance withapplications on comparison of dna and protein sequences.
PloS one , 10(9), 2015.[29] G. TSOs. Netzentwicklungsplan strom 2025. zweiter entwurf der übertragungsnetzbetreiber. , 2015.[30] M. Vega and H. G. Sarmiento. Image processing application maps optimal transmission routes.
IEEE computerapplications in power , 9(2):47–51, 1996.[31] N. West, B. Dwolatzky, and A. Meyer. Terrain based routing of distribution cables.
IEEE computer Applicationsin power , 10(1):42–46, 1997.[32] S. Winter. Modeling costs of turns in route planning.
GeoInformatica , 6(4):345–361, 2002.14
PREPRINT - F
EBRUARY
2, 2021
Appendices
A Accelerating minimal-angle Bellman-Ford for step-functions as angle cost functions
In the proposed minimal-angle Bellman-Ford algorithm (algorithm 1), the edge distances D [ e ] were updated by takingthe minimum edge- plus angle-cost over the incoming edges. At a vertex with k incoming edges and l outgoing edges,the algorithm requires k · l operations to update the distances at all its outgoing edges, since the minimum is taken foreach outgoing edge. Note that these operations appear at each vertex regardless of the type of graph (DAG, directed orundirected). Here, we propose a procedure that decreases the number of update operations significantly if the angle costfunction is a step function: Theorem 2. If c a : [0 , −→ S is a step function with | S | possible angle costs ( | S | < δ − v , δ + v ), the shortest walkcan be computed in time O (cid:0) p · (cid:80) v ( | S | · δ − v + δ + v ) · log( | S | δ − v δ + v ) (cid:1) . As proof, we will explain an algorithm and show its correctness and runtime. To simplify notation, we will from nowon only consider one vertex with k incoming edges and l outgoing edges. We show that in the case of a step-function O ( | S | · k + l · log( | S | kl )) operations suffice to update one vertex (previously k · l ), which leads to the improved runtimeof minimal-angle Bellman-Ford as stated in Theorem 2.We denote the incoming edge vectors by the set E − = { e − , . . . , e − k } and the outgoing edge vectors by E + = { e +1 , . . . , e + l } . Note that if the graph is undirected, it is simply transformed into a directed graph, such that each originaledge appears once in E − and once in E + , and k = l . Thus, from now on a directed graph is assumed. Furthermore,each incoming/outgoing edge reaches/leaves the vertex v in a certain angle. This angle is computed for all edges withrespect to an arbitrary direction (cid:126)w , yielding values α i = (cid:94) ( (cid:126)w, e − i ) for e − i ∈ E − and β j = (cid:94) ( (cid:126)w, e + j ) for e + j ∈ E + .Crucially, we define (cid:94) ( x ) as the angle in clockwise direction from 0 to 360 degrees, instead of a direction-insensitiveangle from 0 to 180 degrees. Only with angle values from 0 to 360 it holds that α i = β j if and only if the vectors e − i and e + j form a straight line ( e + j leaves v in the same direction as e − i reaches it).We further define the angle cost function c a in the following manner: c a ( x ) = θ ≤ x < γ θ γ ≤ x < γ θ γ ≤ x < γ · · · θ | S | γ | S |− ≤ x ≤ In other words, between an angle of γ s − and γ s , the angle cost is θ s . Thus, θ , . . . , θ | S | is the co-domain of the anglecost function c a .In the following algorithm we exploit that if there are only | S | different angle costs, then only k · | S | different values arepossible for the distance of an arbitrary outgoing edge from the source. The set of possible path distances is finite and canbe enumerated as L = { D [ e − ]+ θ , D [ e − ]+ θ , . . . , D [ e − ]+ θ | S | , . . . , . . . , D [ e − k ]+ θ , D [ e − k ]+ θ , . . . , D [ e − k ]+ θ | S | } .The goal is to determine the optimal predecessor for each outgoing edge in E + , by selecting the lowest of these valueswhich is feasible in the respective angle.Our algorithm proceeds in the following manner: First, the angles of the outgoing edges β , . . . , β l are inserted in abalanced and sorted AVL tree ( O ( l log l ) ). Next, the possible path distance values in L are sorted. Sorting | L | = k | S | values requires O ( k | S | · log( k · | S | )) . More specifically, not the distances themselves are sorted, but the correspondingtuples ( i, s ) of the incoming edge index together with the index of the angle cost. The algorithm then iterates over thesorted tuples, such that ( i ∗ , s ∗ ) = argmin i,s D [ e − i ] + θ s is retrieved first. In each iteration, the feasible range of anglesgiven by the bounds α i + γ s − to α i + γ s and α i − γ s to α i − γ s − is considered. All outgoing edges whose anglevalue β lies in this range will have the incoming edge i as its predecessor, if they have not been assigned a predecessorpreviously. 15 PREPRINT - F
EBRUARY
2, 2021In order to retrieve the relevant outgoing angles in this range, we first define an auxiliary function get-inbetween (algorithm 2) that retrieves all values from a tree which lie between a lower and upper bound. Note that this can beefficiently implemented with a recursive function.
Algorithm 2: get-inbetween
Input: AVL tree T , current list of values V , lower bound x l , upper bound x u if T.isEmpty() then return {} n = T.root () // If value was inbetween or too large, go left if x l ≤ n.key () then T l ← sub-tree rooted in left child of n V = V ∪ get-inbetween( T l , V, x l , x r ) // If value was inbetween or too low, go right if n.key () ≤ x r then T r ← sub-tree rooted in right child of n V = V ∪ get-inbetween( T r , V, x l , x r ) // Add tree node if the value is inbetween x l and x r if x l ≤ n.key () ≤ x r then V = V ∪ n return V The main algorithm that was described above is given in algorithm 3.
Algorithm 3:
Step-function update algorithm
Input: E + , E − , α i ∀ i, β j ∀ j // Insert outgoing angle values in AVL tree: T ← empty AVL tree for j = 1 ..l do T.insert( β j ) end // Sort incoming edge distances plus angle costs L = [(1 , , . . . , (1 , | S | ) , . . . , . . . , ( k, , . . . , ( k, | S | )] L ∗ = sort tuples ( i, s ) in L by their distance value D [ e − i ] + θ s // Iterate over list starting with the minimum for ( i, s ) ∈ L ∗ do A ← get-inbetween ( T, {} , α i + γ s − , α i + γ s ) // Add nodes in the other direction A ← get-inbetween ( T, A , α i − γ s , α i − γ s − ) // For all retrieved nodes for β j ∈ A do T.delete( β j ) // update optimal predecessor P [ e + j ] = e − i endend Note that special care has to be taken with respect to the circularity of the x -axis here: As the x -values are angles, thehighest value is actually close to the lowest value in the tree. Thus, if a i + γ s > , then get-inbetween must becalled two times, once with the range of a i + γ s − to 360, and once with 0 to ( a i + γ s ) mod , where mod denotesthe modulo operation. For the sake of simplicity, in algorithm 3 such special cases are not covered, and will not becovered in the following proofs. Similarly, we do not regard the case of several outgoing edges in the same angle (andthus the same β value). Extensions of the above algorithms that regard the border cases and non-unique values arestraightforward to implement. Lemma 7.
Algorithm 3 assigns the correct predecessor to each outgoing edge. PREPRINT - F
EBRUARY
2, 2021
Proof.
Let e + j denote any outgoing edge. If the corresponding angle value β j is retrieved from T , this means β j lies between two values α i + γ s − and α i + γ s due to the correctness of get-inbetween (or between α i − γ s and α i − γ s − , but the proof is analogous). Due to the computation of α and β , the fact that α i + γ s − ≤ β j ≤ α i + γ s implies that γ s − ≤ (cid:94) ( e − i , e + j ) ≤ γ s and the distance thus amounts to D [ e − i ] + c a ( (cid:94) ( e − i , e + j )) = D [ e − i ] + θ s .Furthermore, observe that each outgoing edge is retrieved from the tree only once, since its entry in T is deletedafterwards. When e + j is retrieved from T in an iteration where the current tuple is ( i, s ), then e − i is assigned as itspredecessor and the predecessor is not changed afterwards. For the sake of contradiction, assume there exists a betterpredecessor for e + j , say e − ˆ i . Then this predecessor would imply a lower distance, so D [ e − ˆ i ] + θ ˆ s < D [ e − i ] + θ s withsome other angle value ˆ s . The tuple (ˆ i, ˆ s ) must therefore appear earlier in S ∗ , because the tuples in S ∗ are sorted bytheir corresponding cost. However, if such a better tuple existed, β j would be between the corresponding bounds givenby γ ˆ s − and γ ˆ s , and it would have been deleted from the tree earlier. The contradiction thus shows that the first tuplethat retrieves β j from the tree indeed declares the optimal predecessor for e + j . Lemma 8.
Algorithm 3 requires O (cid:0) ( | S | k + l ) · log( | S | kl ) (cid:1) operations to update the predecessors of all outgoingedges.Proof. The computation of the angle values α , . . . α k and β , . . . , β l is straightforward and requires only one pass overthe incoming and outgoing edges respectively. Inserting all β j ( j = 1 ..l ) in an AVL tree is feasible in O ( l log l ) . Theneach incoming edge is paired with all S possible angle costs and these tuples are sorted, which takes O ( k | S | log k | S | ) steps. In the subsequent loop of k | S | iterations, the nodes are retrieved from the tree, and importantly, once a node wasretrieved it is always deleted. In a single iteration, there could be O ( l log l ) operations, for retrieving and deleting allelements from T . However, the total number of steps, summed over all iterations, is upper bounded by O ( l log l ) aswell because every element can only be encountered once. Since get-inbetween is executed in each iterations (usingat least O (log l ) steps to find the closest values to x l and x r ), there are O ( k | S | log l + l log l ) steps in total in the mainloop.Together, O ( l log l ) + O ( k | S | log k | S | ) + O ( k | S | log l + l log l ) can be summarized as O (( k | S | + l ) log | S | kl ) .This result concludes the proof of Theorem 2, when replacing k by δ − ( v ) , l by δ + ( v ) . Specifically, we use algorithm 1to compute the shortest path, but employ algorithm 3 introduced here as an improved update rule when the anglecost function is a step function. Note that the runtime is only improved if | S | << δ + ( v ) since otherwise | S | δ − ( v ) ≈ δ + ( v ) δ − ( v ) and thus the simple update rule with complexity δ + ( v ) δ − ( v ) would suffice. B An accelerated update algorithm for convex angle cost functions
As in the discrete case above, we show that the runtime decreases for any convex increasing angle cost function. Thealgorithm will be presented in the following.
Theorem 3.
Let c a : [0 , −→ R + be a convex and monotonically increasing angle cost function. Then the runtimeof algorithm 1 reduces to O (cid:16) p · (cid:80) v (cid:0) δ − ( v ) + δ + ( v ) (cid:1) log (cid:0) δ + ( v ) δ − ( v ) (cid:1)(cid:17) , since at each vertex with k incoming and l outgoing edges only O (cid:0) ( k + l ) log kl (cid:1) operations are executed. We define E − , E + and the corresponding in and outgoing angle values α , . . . , α k and β , . . . , β l as above inappendix A. With this notation, we can formulate the problem to determine the the optimal predecessor in terms ofangle and distances for each outgoing edge as ∀ j : Find argmin i D [ e − i ] + c a ( | α i − β j | m ) , | x | m = min {| x | , − | x |} . The operation | x | m computes the smaller angle between both and thereby ensures that the input of c a lies within thefunction domain [0 , .Now D [ e − i ] + c a ( | α i − β j | m ) can be seen as a function associated to the incoming edge i and taking an angle β j as input. We denote these functions by f i such that f i ( x ) = D [ e − i ] + c a ( | α i − x | m ) . Observe that each function isa shifted version of c a ( | x | m ) , displaced by α i in x- and D [ e − i ] in y-direction. Figure 10 shows an example, whereeach incoming edge induces a cost function for the angles of possible outgoing functions. Note that f i is symmetricaround α i , and it has a global minimum at α i because c a is monotonically increasing. Also, if the x-axis is rotated17 PREPRINT - F
EBRUARY
2, 2021 angle value w.r.t. c o s t correspondingpredecessor edge byangle Initial heapFinal heap is predecessor of
Figure 10:
Angle update algorithm for convex, monotonically increasing angle cost functions. The intuition is that the intersectionsof neighboring functions define the borders of the ”area of responsibility”for each incoming edge. Intersections are iteratively addedto a tree, starting with the incoming edge with lowest base cost D . The tree is sorted and balanced. After the tree is fully constructed,each outgoing edge e + j is assigned its optimal predecessor by finding the closest value to β j in the tree. appropriately such that | x | m = | x | then f i is convex because it is a (shifted) composition of two convex functions, c a and | x | . The key observation used in the following algorithm is that - due to convexity - an incoming edge is the optimalpredecessor for a closed range of angles (without gaps). In other words, if two outgoing edges at angles β , β have thesame optimal predecessor, then another outgoing edge with angle β that fulfills β ≤ β ≤ β must have the samepredecessor. To determine the ”area of responsibility”, i.e. the range of angles in which an outgoing edge would get e − i assigned as its predecessor (Figure 10 bottom left), the intersections of the functions f i are computed. The idea is toconstruct a tree T x containing the intersection values (Figure 10 right) that can be traversed in O (log k ) in the end foreach outgoing edge to find the closest intersection and thereby its optimal predecessor edge.First, the incoming edges are sorted by their distance from the source D [ e − i ] . For simplicity of notation, we fromnow on assume that they are already sorted, such that D [ e − ] ≤ D [ e − ] ≤ · · · ≤ D [ e − k ] . Then, the incoming edges aresuccessively processed and the intersection of their respective function f i with previously responsible functions in theirrange is computed and added to a balanced AVL tree T x . algorithm 5 describes the procedure in detail. The tree isbuilt iteratively by insertion and replacement of intersection points. We define an operation compute-intersection that outputs the intersection of two functions. Given f p , f q , compute-intersection ( f p , f q ) yields a x value at whichthe functions intersect. Crucially, the angle values are circular, i.e. an angle of 0 degrees corresponds to an angleof 360 degrees. Therefore, compute-intersection ( f p , f q ) (cid:54) = compute-intersection ( f q , f p ) because if the twofunctions intersect at any point x , then they intersect at least at two points, x and − x . Thus, for each intersectionwe store which function was lower in clockwise (positive x-) direction and which one in counter-clockwise (negativex-) direction. For example, a triple ( x, p, q ) means that e − p is the optimal predecessor for β -angles that are close to x in counter-clockwise direction and e − q is optimal for outgoing edges close to x in clockwise direction. Due to thecircularity it is for example possible that α p = 300 , α q = 10 and x = 330 . It will be argued further below that compute-intersection can be approximated efficiently.Secondly, algorithm 5 utilizes an operation called find-closest that returns the two closest value in a balancedAVL-tree. The algorithm is outlined in detail in algorithm 4. The balanced and sorted tree is traversed and the closestleft and right values σ l , σ r with respect to the input element are memorized and returned.Last, in addition to the sorted tree of intersections T x we also maintain a double-linked list L of neighboring intersections.If ( x, p, q ) is the intersection triple of f p and f q , then ( x, p, q ) .next is defined as the next intersection in clockwisedirection which must therefore be of the form ( y, q, r ) with some r . Similarly, ( x, p, q ) .previous is of the form ( x, o, p ) and denotes the next intersection counter-clockwise. In the main algorithm detailed in algorithm 5, foreach new incoming edge e − i the closest angle values of incoming edges that were already processed are found with18 PREPRINT - F
EBRUARY
2, 2021 find-closest . To do this efficiently in logarithmic time, all α i are iteratively inserted into a second tree T α and find-closest ( T α , α i ) is called. Algorithm 4: find-closest
Input: AVL tree T element x // define σ l and σ r as the currently closest elements σ l , σ r = ∞ n = T.root () // Traverse tree while not isLeaf ( n ) doif x > n.key () then σ l = nn = getRightChild ( n ) else σ r = nn = getLeftChild ( n ) end // check border cases if one is still infinity if σ l is ∞ then // set left bound to the largest element in the tree σ l = T [ − if σ r is ∞ then // set right bound to the smallest element in the tree σ r = T [0] return σ l , σ r In short, algorithm 5 implements the following steps: At the i -th step, f i is retrieved, thus corresponding to the i -thlowest incoming cost D [ e − i ] . Then find-closest is called for a i . Let the result, i.e. the closest angles to α i , bedenoted by α p and α q . If f i obtains a better cost at the intersections between f p and f q , the triple is deleted from T x and the next triples in both directions are checked. This is repeated until encountering an intersection point where f i can not improve on the previous functions. Due to convexity, once f i is outperformed by another function it will not belower again in this direction.Once these functions that mark the end of the area of responsibility of edge i have been determined, the intersectionswith them are computed and added to T x . Also, the double-linked list is updated accordingly.After all incoming edges were considered and possibly added to T x , the outgoing edges can be processed. For each β j ,the left and right closest intersection triples ( x , − , i ) and ( x , i, − ) are found with a call to find-closest and thecorresponding incoming edge e − i is assigned as the predecessor of e + j .Consider the example in Figure 10. First of all, the first function must be found that intersects with f (PART 1 inalgorithm 5). Clearly, f ( α ) ≤ f i ( α i ) ∀ i , so it must be checked whether f is better at any other point, formallywhether ∃ x : f ( x ) > f ( x ) . Lemma 9 will show that it is sufficient to check whether f is lower at the opposite point of α , namely ( α + 180) mod , where f is maximal. Since f satisfies this condition, the two intersection (clockwiseand counter-clockwise from α ) between f and f are computed and added to the tree as shown in Figure 10 (upperright). T x is sorted and balanced. In the first iteration of PART 2 in algorithm 5, α is considered and its neighboring α i are determined with find-closest . In the example, α is between α counter-clockwise and α clockwise. Wetherefore check whether f improves the value at intersection x , . We indeed find that f ( x , ) < f ( x , ) , so thereis clearly a range of angles where e − is the optimal predecessor. The triple ( x , , , is thus deleted from T x . Todetermine the range where f is optimal, the algorithm 5 further proceeds by checking the next intersections in clockwiseand counter-clockwise directions. In this case, the next intersection in both cases is x , . However, f ( x , ) > f ( x , ) ,so the algorithm can stop here. After these checks, in algorithm 5 now p ∗ = 2 for our example, and q ∗ = 1 . We thuscompute the new intersections x , and x , and add the triples ( x , , , and ( x , , , to T x . On the other hand, inthe last iteration f (whose angle α lies between α and α ) does not outperform the value at the intersection between f and f and consequently α is not an optimal predecessor for any outgoing angle. After all incoming edges were19 PREPRINT - F
EBRUARY
2, 2021considered and T x is final, the β -values of the outgoing edges are processed. In Figure 10 for example β is betweenthe intersection triples ( x , , , and ( x , , , and thus e − is the optimal predecessor for e +3 . Algorithm 5:
Efficient angle cost update algorithm for convex angle cost functions // PART 1: Find first intersecting functions i start = 2ˆ α ←− ( α + 180) mod // opposite point to α while i start ≤ k and f i start ( ˆ α ) > f ( ˆ α ) do i start ←− i start + 1 endif i start > k then // Break: no function better than f at any point T x .insert( α , , ) return T x x = compute-intersection ( f , f i start ) x = compute-intersection ( f i start , f ) ( x , , i start ) .next = ( x , , i start ) .previous = ( x , i start , x , i start , .next = ( x , i start , .previous = ( x , , i start ) T x .insert( ( x , , i start ) , ( x , i start , ) T α .insert ( α , α i start ) // initialize T α // PART 2: Process further incoming edges for i = i start ..k do // Get closest α values α p , α q ←− find-closest ( T α , α i ) x p,q ←− compute-intersection ( f p , f q ) // Does f i outperform the previous functions at one or more intersections? if f i ( x p,q ) < f p ( x p,q ) or f i ( x p,q ) < f q ( x p,q ) then T α .insert ( α i ) T x .delete( ( x p,q , p, q ) ) // Check counter-clockwise intersection ( x , o, p ) ←− ( x p,q , p, q ) .previous while f i ( x ) < f p ( x ) do // delete and get next counter-clockwise T x .delete( ( x , o, p ) ) T α .delete( α p ) ( x , o, p ) ←− ( x , o, p ) . previous end p ∗ ←− p // Check clockwise intersection ( x , q, r ) ←− ( x p,q , p, q ) .next while f i ( x ) < f q ( x ) do // delete and get next clockwise T x .delete( ( x , q, r ) ) T α .delete( α q ) ( x , q, r ) ←− ( x , q, r ) .next end q ∗ ←− qx ∗ ←− compute-intersection ( f p ∗ , f i ) x ∗ ←− compute-intersection ( f i , f q ∗ ) T x .insert( ( x ∗ , p ∗ , i ) , ( x ∗ , i, q ∗ ) ) ( x , o, p ∗ ) .next = ( x ∗ , i, q ∗ ) .previous = ( x ∗ , p ∗ , i )( x ∗ , p ∗ , i ) .next = ( x , q ∗ , r ) .previous = ( x ∗ , i, q ∗ ) end // PART 3: process outgoing edges for j = 1 ..l do ( x , p, q ) , ( x , q, r ) ←− find-closest ( T x , β j ) // update predecessor P [ e + j ] = e − q end In the following, we will describe several properties of the proposed method that will be necessary to proof thecorrectness and runtime efficiency of algorithm 5. 20
PREPRINT - F
EBRUARY
2, 2021
Lemma 9.
Let f p and f q be functions for the incoming edges e − p and e − q as defined above, and let D [ e − p ] ≤ D [ e − q ] and α p < and α p < α q < α p + 180 , such that f q is equal to f p shifted in positive x -direction and in positive y -direction. Let ˆ α p = α p + 180 be the point on the opposite side of α p on the cycle.If ∃ x : f q ( x ) < f p ( x ) then ∀ x (cid:48) ∈ [ x, ˆ α p ] : f q ( x (cid:48) ) < f p ( x (cid:48) ) Proof.
We will use a property of convex monotonically increasing functions that can be seen as an increasing returnsproperty:
Claim 1.
For a convex function g and a > b , a, b, c ≥ it holds that g ( a + c ) − g ( a ) ≥ g ( b + c ) − g ( b ) (B.1) Proof.
The intuition is that if g is differentiable its gradients are increasing and thus the difference in value becomeslarger. To prove the claim, we use another well-known property of convex functions, namely if f is convex and f (0) ≤ then ∀ d, c ≥ f ( c ) + f ( d ) ≤ f ( c + d ) (B.2)We define f ( x ) = g ( x + b ) − g ( b ) . Since f is only a shifted version of g , f is convex as well and also note that f (0) = g ( b ) − g ( b ) = 0 . We further set d = a − b > and plug in d in property B.2 to yield f ( d ) + f ( c ) ≤ f ( d + c ) ⇐⇒ f ( a − b ) + f ( c ) ≤ f ( a − b + c ) ⇐⇒ g ( a − b + b ) − g ( b ) + g ( c + b ) − g ( b ) ≤ g ( a − b + c + b ) − g ( b ) ⇐⇒ g ( a ) − g ( b ) + g ( b + c ) ≤ g ( a + c ) which is equivalent to the claim.To continue the proof of lemma 9 we consider the point x where f q ( x ) < f p ( x ) . Observe that due to D [ e − p ] ≤ D [ e − q ] and with c a monotonically increasing, f q ( x ) < f p ( x )= ⇒ D [ e − q ] + c a ( | x − α q | m ) < D [ e − p ] + c a ( | α p − x | m )= ⇒ | x − α q | m < | α p − x | m Let x (cid:48) be any other point between x and the opposite point of α p as defined above ( x (cid:48) ∈ [ x, ˆ α p ] ). Since c a ismonotonically increasing, f p is maximal for ˆ α p . We use the claim show that f q ( x (cid:48) ) < f p ( x (cid:48) ) by setting a = | x − α p | m and b = | x − α q | m (leading to a > b ), and further setting c = | x − x (cid:48) | m : f q ( x (cid:48) )= D [ e − q ] + c a ( | x (cid:48) − α q | m ) ≤ ∗ D [ e − q ] + c a ( | x (cid:48) − x | m + | x − α q | m )= D [ e − q ] + c a ( c + b ) ≤ B. D [ e − q ] + c a ( c + a ) − c a ( a ) + c a ( b )= D [ e − q ] + c a (cid:0) | x (cid:48) − x | m + | x − α p | m (cid:1) − c a (cid:0) | x − α p | m (cid:1) + c a (cid:0) | x − α q | m (cid:1) = ∗∗ D [ e − q ] + c a (cid:0) | x (cid:48) − α p | m (cid:1) − c a (cid:0) | x − α p | m (cid:1) + c a (cid:0) | x − α q | m (cid:1) = f q ( x ) + c a (cid:0) | x (cid:48) − α p | m (cid:1) − c a (cid:0) | x − α p | m (cid:1) < f q ( x ) EBRUARY 2, 2021(*) follows from the triangle inequality that also holds in the circular case. (**) is based on the condition that x (cid:48) ∈ [ x, ˆ α p ] .In short, due to the property of increasing returns, if f q is lower than f p at any point it will remain lower until intersectingwith f p from the other side and thus at least until ˆ α p .In Lemma 9 it was assumed that α p < α q < α p + 180 which might seem very restrictive. The assumption was madeto reduce complexity in the proof, but the statement also holds for the circular case. To see this, observe that giventwo functions f p and f q , the x -axis can be rotated and flipped to fulfill the above condition. Since it is circular, thisoperation is always possible and does not change the shown property. The same holds for the following observation,that formally shows that the any new function must only be evaluated at the intersection point of the previous functionsin order to determine whether it is superior anywhere. Lemma 10. Let α i refer to an angle of an incoming edge between α p and α q , such that in the simplified non-circularmodel it can be assumed that α p < α i < α q . Assume also that D [ e − i ] ≥ D [ e − p ] and D [ e − i ] ≥ D [ e − q ] . Let x p,q denotethe intersection of f p and f q .If f i ( x p,q ) > f p ( x p,q ) = f q ( x p,q ) then (cid:64) x : f i ( x ) < f p ( x ) ∧ f i ( x ) < f q ( x ) .Proof. We again use an observation that was made in Lemma 9: If D [ e − i ] ≥ D [ e − p ] , then f i ( x ) < f p ( x ) implies | x − α i | m < | α p − x | m and the same holds for q . Thus, if there exists an x where f i ( x ) < f p ( x ) and f i ( x ) < f q ( x ) then it follows that α p < x < α q . Assume α p < x ≤ α i < α q (the other case is analogous). Due to Lemma 9if f i ( x ) < f p ( x ) at x then it remains lower until the opposite point ˆ α p . Since the intersection x p,q is at ˆ α p at thelatest (otherwise f q would not be better then f p anywhere), it can be concluded that f i outperforms f p also at theintersection. Theorem 4. algorithm 5 assigns the correct predecessor to each outgoing edge.Proof. Let e − i be the optimal predecessor for an outgoing edge j ∗ . In the first loop of algorithm 5 an incoming edge e − i is only skipped if its value at ˆ α is higher than for e . By Lemma 9, if the value is not better than f at ˆ α , then it islarger everywhere. Consequently, any incoming edge that is skipped in the first loop of algorithm 5 is outperformed by e − and would thus not be the unique optimal predecessor of any outgoing edge.Any e − i that is not skipped in the PART 1 will be considered at some point in PART 2 in algorithm 5. Its closest angles α p and α q are retrieved. Due to Lemma 10 it suffices to check whether f i is better at the intersection of f p and f q ,otherwise it can be discarded. Lemma 10 can be applied because if α p and α q were retrieved from T α , then they havealready been processed earlier than α i and consequently their distance D must be lower.Furthermore, if f i is indeed lower at the intersection, then there might be other intersections where f i outperformsprevious incoming edges. Note that any of such points must lie between α p and α q . algorithm 5 takes care of suchspecial cases by checking the neighboring intersections in both directions. If at some point it is outperformed by anotherfunction, then it can not become better at any further point due to convexity. This procedure thus defines the currentarea of responsibility for edge e − i .With further iterations, this area can be reduced in size or removed. Due to the conditions in algorithm 5, it is onlyremoved if a new incoming edge indexed with ˆ i obtains lower values at both intersections of f i with its neighboringfunctions. Then algorithm 5 would correctly compute the new intersections of f ˆ i with f p and f q which replace theintersections of f i . It would thus not be possible that f i is the optimal predecessor of any outgoing edge, because it waspreviously outperformed behind the borders of its area of responsibility by the respective intersecting functions, andnow it is also outperformed in the area inbetween.In the last for-loop (PART 3) in algorithm 5 concerning the outgoing edges, the closest intersections ( x , u, v ) and ( x , v, w ) to β j ∗ are retrieved from T x . Note that by the procedure of algorithm 5 it is ensured that find-closest always yields triples of this form where v appears in both triples and ( x , u, v ) . next = ( x , v, w ) and ( x , v, w ) . previous = ( x , u, v ) . To see this, observe that the while-loops in algorithm 5 stops at two triples ( x , q ∗ , r ) and ( x , o, p ∗ ) . Since all intersections inbetween were deleted, it is correct that ( x , q ∗ , r ) is assigned thenew intersection with of f i with f q ∗ as the previous triple and equivalently ( x , o, p ∗ ) is assigned the new intersectionof f i with f p ∗ . Lemma 11. The compute-intersection operation can be approximated efficiently and sufficiently for algorithm 5. PREPRINT - F EBRUARY 2, 2021 Proof. While the efficiency to compute an intersection between a function and its shifted version depends on thespecific function, it is actually not necessary here to know the exact intersection. Crucially, in the end we are onlyinterested in the values β , . . . , β l , so it suffices to evaluate the functions at these values. Given two functions f p , f q andtheir theoretical intersection x p,q , we just aim to find j ∗ = argmin j | x p,q − β j | . For this purpose a similar procedure asthe find-closest function (algorithm 4) is utilized. First, all β j ( j ∈ [1 ..l ]) are inserted in an AVL tree T ( O (cid:0) l log l (cid:1) ).For each call of compute-intersection , the tree is traversed in the following manner: Starting at the root node, f p and f q are evaluated at the β value stored in the node. If f p ( β ) > f q ( β ) , the left child node is selected, otherwise theright child node. The intuition is that if f p ( β ) > f q ( β ) then the intersection must be at a smaller β value where f p takeslower values. During traversal of the tree, the optimal node with respect to the distance | f p ( β ) − f q ( β ) | and its distanceitself are memorized such that the closest β is determined once a leaf node is reached.However, one might object that the border case of the circularity is problematic in this case. Indeed, it might happen thatthe rightmost leaf in the tree is reached but the optimal β value (the one closest to the intersection) is not the leftmostone. To prevent this case, we propose to build two trees with all β j inserted: One sorted and balanced tree with all β j ,and one with all ˆ β j where ˆ β j = ( β j + 180) mod . After both trees have been traversed, the value that is closer tothe intersection, i.e. where | f p ( β ) − f q ( β ) | is lower, is selected. This ensures that the closest β value is found, as it isnot a border case in at least one of the trees.We claim further that the discrete evaluations do not impact the proof of Theorem 4. Assume that α p ≤ β p,q ≤ x p,q (the other side is again analogous). Note that there can not exist any ˆ β that satisfies α p ≤ β p,q ≤ ˆ β ≤ x p,q becausethen ˆ β would have been selected as the intersection between f p and f q .On the other hand, there can be β (cid:48) with α p ≤ β p,q ≤ x p,q ≤ β (cid:48) ≤ α q . Crucially, it is possible that f i ( β (cid:48) ) < f p ( β (cid:48) ) and even f i ( x p,q ) < f q ( x p,q ) , but f i ( β p,q ) > f q ( β p,q ) such that the new function f i is not better at the approximatedintersection but still lower at the actual intersection. However, it has been shown in Lemma 9 that in this case f i ( x ) < f q ( x ) ∀ x < β (cid:48) . This is why in algorithm 5 it is checked whether both f p and f q outperform f i at theapproximate intersection (which would be redundant if the exact intersection was available because the values for f p and f q would be equal). Thus, if β (cid:48) with the above properties exists, then f i ( β p,q ) < f q ( β p,q ) and therefore thealgorithm would still compute the intersections of f i with the previous functions and β (cid:48) would be assigned to the correctrange. Theorem 5. The runtime of algorithm 5 is O (cid:0) ( k + l ) log kl (cid:1) .Proof. First, the incoming edges are sorted by their distance to the source vertex, taking O (cid:0) k log k (cid:1) time. As argued inlemma 11, the compute-intersection operation can be approximated here by a discrete evaluation on the β -values,leading to a runtime of O (cid:0) log l (cid:1) for each call to compute-intersection once the AVL-tree of β -values is build( O (cid:0) l log l (cid:1) ). In addition, each insert and delete operation to T x takes O (cid:0) log k (cid:1) since at most k elements are in thetree (in each iteration, one intersection is deleted and two intersections are added). Similarly, each insert and deleteoperation, as well as each call to find-closest on T α requires O (cid:0) log k (cid:1) time. Last, the find-closest operationagain only requires a traversal of T x in O (cid:0) log k (cid:1) steps.Now it is important to observe that in the while-loop of PART 2 of algorithm 5, elements are only deleted from T x and T a . The number of operations in the while loop is therefore limited to the maximum number of elements of the tree.The tree contains at most O ( k ) values at any point, because in each iteration (for each incoming edge) at most two oneintersections are added, and one is deleted. Thus, PART 2 of algorithm 5 requires O (cid:0) k (log l + log k ) (cid:1) = O (cid:0) k log kl (cid:1) operations.Finally, in PART 3 of algorithm 5 the closest element in T x is found for each outgoing edge, causing an additionalruntime of O (cid:0) l log k (cid:1) . Given the optimal predecessor edges, the distances of E + can be updated in O (cid:0) l (cid:1) . Altogether,the runtime is O (cid:0) k log k + l log l + k log kl + l log k (cid:1) ∈ O (cid:0) ( k + l ) log kl (cid:1)(cid:1)