Computing the trajectories for the development of optimal routes

Planning the construction of new transport routes or power lines on terrain is usually carried out manually by engineers, with no guarantee of optimality. We introduce a new approach for the computation of an optimal trajectory for the construction of new transit routes and power lines between two locations on a submanifold U⊂R3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U\subset \mathbb {R}^{3}$$\end{document} representing the topography of a terrain. U is approximatively modeled by a special weighted grid. On this grid, the shortest paths for the construction of new routes are determined, whereby we consider three optimization criteria: routes with minimum distance, routes with lowest construction costs and routes with minimum absolute altitude variations or minimum absolute gradients. Subsequently, a combination of these criteria is used to expand this problem into a multi-criteria optimization problem. A shortest path algorithm, such as the Dijkstra algorithm, is used to compute optimal compromises for the construction of new routes.


Introduction
Worldwide, many developing countries have an undeveloped technical and logistical infrastructure. This not only causes numerous problems for the inhabitants in their daily life, but also has negative economic consequences (more on this can be found in Fay and Yepes (2003), Gaal (2017), Gurara et al. (2017), Oyedele (2016)).

3
Classical examples of a logistical infrastructure are transport networks such as roads, tracks, and waterways. On the technical side, these are supply networks such as power lines, freshwater and wastewater lines, and communication networks. An immature technical infrastructure can severely curb the energy and transport sector. These containments can not only increase the manufacturing and transport costs of a product but can also impact entire construction projects, transport networks, and economic networks. The resulting consequences can reduce the development and progress of these countries. Besides, investors lose interest in investing in underdeveloped countries with a high amount of natural resources since the risk of a misinvestment is too high. Infrastructure projects are usually very costly. Planning firms or consulting companies often plan the development of new routes (alignment) based on empirical values or using estimation methods. However, these methods do not guarantee optimality. Such planning projects and feasibility studies can not only cause high costs but also be very time-consuming. In the case of extensive studies in risky areas, the presence of the staff on-site can pose a threat to the safety of the personnel. Risky zones can include regions with natural disasters or conflict zones. Last but not least, manual planning involves a residual risk of corruption and manipulation, as a manual planned route allows for unnecessary detours. By planning and constructing redundant detours, companies can generate more turnover and use part of the turnover to influence the transport ministries or the relevant authorities.
The Shortest Path Problem is a well-known graph theory problem where the goal is to find a shortest path between two nodes (start and end points) within a weighted graph G = (V, E, w) , where V is a set of nodes, E is a set of edges, and w ∶ E ⟶ ℝ + is a cost function that defines the weight w i,j of an edge {i, j} ∈ E . A shortest path between two nodes is a route with a minimum sum of all used edge weights to get from the start node to the destination node. Shortest path algorithms are particularly used for path planning or finding a path. Route planning is often used for autonomous mobile robots to find an optimal route from a starting point to a destination in complex environments. There are different heuristic and non-heuristic approaches for finding a way to calculate an ideal path on a terrain. One approach is that a graph is given and a search algorithm for a path is performed on this graph. Another approach is that a certain search algorithm does not search within a given graph, but simultaneously creates a graph in the direction of the target on which a potential optimal path could lie while searching for a path. Mitchell and Sharir present in Mitchell and Sharir (2004) an efficient algorithm that calculates Euclidean shortest paths in three dimensions with obstacles, in polynomial time. An L 1 shortest path on polyhedron terrain, L 2 shortest path on parallel walls, and L 2 shortest path on a landlike, layered, and axis-aligned rectangles are analyzed and computed. Belaidi et al. describe in Belaidi et al. (2014) an algorithm for mobile robots in a polygonal grid that plans a path to cross a given terrain with irregularities and obstacles. The objective function is, on the one hand, the distance and on the other hand the speed crossing an area. Saranya et al. modify the D * algorithm in Saranya et al. (2016) by defining a cost function that depends on the distance and slope of the terrain. This allows an autonomous mobile robot to calculate an optimal route between two points in an unknown environment, such as an irregular terrain. In Muñoz et al. (2017), Muñoz et al. present another algorithm, the so-called "3D Accurate Navigation Algorithm (3Dana)", which finds the most convenient path for mobile robots on a terrain (DTM-Digital Terrain Model) by considering, for example, the maximum gradient of a terrain that can be driven on. Subsequently, they tested this algorithm on the digital topography of Mars. The geometry of a terrain influences the trajectory of a potential route. Fowler and Little show in Fowler and Little (1979) that they model a terrain with a Triangulated Irregular Network (TIN). In this work they describe a method to automatically extract a TIN model from dense raster data. A first approximation is constructed by automatically triangulating a set of feature points derived from the grid model. Basu and Snoeyink describe in Basu and Snoeyink (2007) how RTINs (rectangular, irregular networks) can be used to represent terrain as a mosaic of polynomial surface areas.
Compared to the above papers, the main focus of our research is not the development of an algorithm for finding a shortest path, but the development of algorithms for a discrete modeling of a Riemannian submanifold U ⊂ ℝ 3 through a special grid = (V, E, W) of nodes V, edges E and weights W, where U is a topographic terrain with obstacles. On a shortest path for the expansion of transport infrastructure on U is searched.
A regular grid is constructed on a geographical area A by mapping to the topography U of A and representing U approximatively. A family of grids a , b , c with different weights W a , W b , W c characterized by certain criteria is calculated. Each edge {i, j} represents a segment of a route, and the non-negative weight w i,j of that edge indicates (a) the length, (b) the construction cost, or (c) the absolute elevation change between nodes i and j. For the calculation and evaluation of the edge weights W a , W b , W c different algorithms were developed. For these algorithms different data are used, e.g., the soil condition, the elevation data of a terrain (topography) and the construction costs per unit length. The resulting routes can be either transport routes or transmission lines connecting two locations or running from one starting point through several locations to a destination. The aim is to construct the shortest route in terms of distance, the most costeffective route, and the most convenient route with minimum vertical fluctuation. Compared to the above approaches, we will not use heuristics to search for an optimal route but solve the problem with a search algorithm, such as the Dijkstra (1959), which guarantees a global optimum. Through the Pareto method, these different grids are combined into multi-criteria grids Pareto = Pareto 1 , … , Pareto n so that a family of n grids is generated based on compromises between these criteria (Ehrgott 2005). On each multi-criteria grid Pareto k of Pareto an optimal route is determined. These solutions represent a set of optimal trade-offs to design ideal transport and energy infrastructures. The presented work is based on Zazai (2020) and is a largely extended version of Zazai and Fügenschuh (2018).
The advantages and differences of our method compared to the above approaches in the papers are as follows: 1. We present a new method by planning infrastructure development with discrete optimization and graph theory, considering different construction criteria. This approach is different from the approach of autonomous mobile robots, where a vehicle moves from its starting point to its destination in an environment. The environment remains invariant while an autonomous mobile robot crosses it, but when a route is laid out, the environment may change as a result of construction work. This means that certain topographical outliers can be tolerated because they are removed when a route is developed. 2. Although the TIN (Triangulation) method is a good approximation for modeling a terrain, it does not guarantee efficiency for calculating routes on that terrain. Especially when the geometry of a terrain is only a secondary consideration, such as when designing the most cost-effective route, where the ground conditions must be taken into account. 3. Unlike the path-finding of autonomous mobile robots, the edges of a grid can consist of different infrastructure segments, such as roads, tracks, pipelines or power lines. For the selection of a specific infrastructure segment, the input data is modified and adjusted, such as the average cost or the structure of the ground. 4. On no locally, but globally optimal routes are computed. 5. Once a grid of a topographic submanifold U is computed, it is possible to calculate between any start and destination point on U an optimal path with a running time of O(|V|log(|V|) + |E|)Turau (2004). For each changed start and destination point on U, the cost for the computation of the edge weights W is being saved. 6. All planning results, such as trajectory coordinates, construction costs, distances, and total absolute altitude changes, are given separately for each optimal route. 7. The Pareto method is used to define multi-criteria grids by using the different grids a , b , and c . The results can model an environment more realistically. These grids are used to determine various optimal trajectories with the corresponding results, which represent optimal compromises between the mentioned criteria.

Topographical grid generation
In this section, we construct a special grid on a spheroid S. For this we take an area U ′ from S, which is curved. This curved area has to be mapped on a Euclidean plane in ℝ 3 by a homeomorphism. Thus the construction of in the Euclidean space will be simplified. A spheroid S can be described by an implicit representation or a parameter representation where the numbers a and c are the half-axes of the rotating half-ellipse, − 2 ≤ ≤ 2 , 0 ≤ ≤ 2 and a, c > 0.
The shape of the earth is spheroidal which defines an equivalent manifold to S. In the following, S represents the surface of the earth. Since the earth is represented by a flattened spheroid, we set a > c . A submanifold U ′ of S is mapped to a submanifold U of ℝ 3 by a homeomorphism U ′ so that U is a Euclidean representation of the curved surface U ′ in ℝ 3 . It is assumed that U describes the topography of a geographical area A ⊂ (x, y, 0) ∈ ℝ 3 on Earth. So the surface of U can have irregular altitudes and depths. The topography of an area is collected from discrete elevation data which is given by a matrix H. A rectangular weighted grid of nodes and edges is constructed on a geographical area A, so that (A) defines a geographical grid of A. (A) consists of rectangular grid cells that have vertices and borderlines in the plane in ℝ 3 . The boundary of each cell consists of four borderlines. The upper and lower borderlines consist of m ∈ ℕ equidistant edges and the left and right borderlines consist of k ∈ ℕ equidistant edges (see Fig. 1). A cell of has 2m + 2k associated nodes on the borderlines. All node pairs of these cells (except those on the same borderline) are connected by edges. These edges are added to the set E A . A grid consists of X ⋅ Y cells, where Y cells in are stacked vertically and in every i− th row X cells are placed horizontally next to each other.
A grid with rectangular cells has the advantage that it is a structured grid. The cells are regularly sorted (parqueted) so that they lie next to each other in a grid. This allows the nodes of the cells to be easily indexed and counted one after the other. All node pairs of a rectangular cell of with 2k vertical and 2m horizontal segments are connected by m 2 + 4km + k 2 edges. On the surface of U, a grid (U) = V U , E U is constructed from nodes V U and edges E U , which approximatively models the topography of A. Next, the topographic grid is to be extended with graphs from altitude contours that lie on local areas of submanifold U, such as the contour line of a mountain or the altitude edges of a canyon. These contours improve the modeling quality of terrain. Since the elevation data H are discrete, only discrete contours can be generated from these data.
Calculating these contours for each altitude h i ∈ h = h 1 , … , h k of U the Marching Squares algorithm (see Ho et al. 2005;Hanisch 2004) is used, where k is the number of contours for different altitudes from H. Consequently, the result is a set of edge lines in different altitudes of h. Each contour consists of nodes and Each contour defines a path of edges on U. The contours are required later for the computation of a mountain road. These contours which are adjacent and lie on top of each other are interlinked by additional edges, so that one has a transition between the contours in different levels. As an illustration, two-layered contours with level h i and h i+1 are considered which are connected as follows: Let h i = V i , E i and h i+1 = V i+1 , E i+1 be two adjacent contours lying one above the other that generate a graph through nodes V i , V i+1 and edges j is a point representing a node u j ∈ V u and dist is a distance function between any two points in ℝ 3 . Each node v j is connected to its nearest neighboring node u j through an additional edge e j = v j , u j . Repeat this step for all contours except the second last. This generates a mesh between neighboring contours that defines an unstructured grid. Figure 2 shows two neighboring contours h i and h i+1 where the red edges connect all nodes of V i to the nodes of V i+1 . Some nodes of V i+1 located far away have no direct connections to the nodes of V i . The crosshatched red edges make a connection to the next nodes of the neighboring contours.
In the last step, these contours are connected to the grid (U) = V U , E U so that a transition between them is possible. To connect (U) with the set of contours , k} , edges with minimal distance are searched, which connect these two sets. As far as possible, all nodes of the contours should be connected with V U . Seven cases are considered for this type of meshing to keep the number of connections to a minimum except for the most necessary: Let h i−1 , h i and h i+1 be three contours lying on top of each other.
(  Figure 3 shows as an example three discrete height contours of U, which on the one hand are connected by dashed red edges and on the other hand these contours are connected to the (red) nodes of the grid (U) = (V U , E U ) by orange edges. The dashed blue connections are the edges that connect the nodes of V U . The minimum distance from v j ∈ V i to V U can be determined for the above seven cases as follows: Let p i j be an associated point of v j ∈ V i , and (4) two subsets of V U and q 1 , … q |V + i ∪V − i | points in ℝ 3 representing nodes of these subsets. To minimize the distance of each point p i j of the node v j ∈ V i the following applies with the auxiliary condition Repeat this algorithm for each node of the contours of U until they have meshed with the grid (U) . This mesh, together with the mesh of all neighboring and superimposed contours of the set (h) , generates an unstructured mesh ( ) = V , E consisting of nodes V and edges E . Since (U) is structured, together with ( ) it forms a hybrid grid This grid is used in the next section to calculate optimal routes or shortest paths to the surface of U.

Optimal routes
The construction of an optimal route on a topographical area U is based on the computation of a shortest path. To solve this problem, the routes are discretized to a grid . The edge weights of are non-negative values, which are determined by a cost function w ∶ E ⟶ ℝ ≥0 , e i,j ⟼ w i,j . Since these values are always non-negative, algorithms such as the Dijkstra (1959) can be used to compute a global optimal route that represents the shortest, most cost-effective, or most convenient route (minimum absolute elevation changes along the route). The higher the number of cells is chosen at a constant size of , the more accurately the routes are determined with simultaneously higher computing time.
Obstacles on U can be modeled by finitely many discrete polygons O 1 , … , O m . If a trajectory is not allowed to cross an obstacle, then it must first be checked whether an edge e i,j intersects one of the polygons O 1 , … , O m . If e i,j has at least one intersection with these polygons, then the weight of this edge is set as w i,j = +∞ . If e i,j touches these polygons at the boundary or has no intersection point, then the weight of this edge is given by the equation (12). If a set of obstacles are O ∶= O u , … , Q k for u ≥ 1 and k ≤ m lies between a starting node and a destination node then a shortest route is computed around O u , … , Q k .
The geographical coordinates on earth are usually given in latitude and longitude coordinates. The latitude coordinate is the geographical width measured from the equator to the north and south. The latitude values range from −90 • at the south pole over 0 • at the equator to +90 • at the north pole. The longitude coordinate is the geographical length going from the zero meridian 0 • to 180 • eastwards and from 0 • to −180 • westwards.

Shortest routes
Let S be a spheroid and U S a submanifold of S, where S has the size of the earth and the subset U S has a regular non-topographic surface. A submanifold U is defined as a connected topographical area on earth which represents the submanifold U S in ℝ 3 as a topographical surface. The shortest routes with respect to a distance are routes within a weighted topographic grid ℝ 3 (U) = V, E, W D , where the weights W D indicate the distances of the edges of E U . First W D is computed and then the weights of W D are assigned to the corresponding edges of the grid ℝ 3 (U) . For U S and U the latitude and longitude coordinate systems are used. The distances W D of the edges of ℝ 3 (U) are computed with the formulas of Vincenty for spheroids (Vincenty 1975) since U is an approximation of U S except for the topography. Let be a function that gives a distance between two points p i and p j of U according to Vincenty formulas for distances on U S . Due to the topological conditions, the edges of a topographic grid ℝ 3 (U) have different slopes. Usually, traffic routes may not exceed a certain slope otherwise, the crossing of this route is not possible or the consumption of a vehicle increases significantly. Therefore, when building a route, a maximum gradient angle is always taken into account. An edge where v i,j (k) for k = 1 … u i,j − 1 are elements of a set of auxiliary nodes V aux lying on the edge e i,j (see Fig. 4). For the first and last auxiliary node of an edge e i,j is gives the altitude of v i,j on an area A ⊂ (x, y, 0) ∈ ℝ 3 where A has the altitude 0. The absolute elevation change between two neighboring auxiliary nodes v i,j (k − 1) and v i,j (k) of V ∪ V aux is defined as for k = 1, … , u i,j . A segment e i,j (k) is called critical if this segment has at least one intersection with O so that this segment is either entirely or partially in O. If e i,j (k) is a critical segment then the edge e i,j is also critical. All critical segments on ℝ 3 (U) with respect to O define a set E crit . Each segment e i,j (k) is assigned a distance d i,j (k) : Let ∶ U ⟶ A be an orthogonal projection of the typographic area U on the area A. With Vincenty's formulas,first the length of the segment e i,j (k) is calculated where the points p i,j (k − 1) and p i,j (k) represent the nodes v i,j (k − 1) and v i,j (k) for k = 1, … , u i,j in ℝ 3 . Finally, the length of a segment is computed by considering a absolute slope angle or maximum slope angle that can cause an additional lengthening: for all k ∈ 1, … , u i,j , max > 0 and N ∈ ℝ + . Serpentines are mostly built in mountain slopes to be able to go up or down steep areas. The second entry in equation 12 describes the segment e i,j (k) as a serpentine with an maximal slope angle max . That is, if this segment exceeds the angle max , then this segment is assigned the length of a serpentine which has a maximum slope angle max . A serpentine is due to its shape longer than a straight segment. If a segment e i,j (k) is within one of the polygons of O or it intersects with O then the value N in the third entry of equation 12 is the weight of this segment. The polygons of O are usually unfavorable areas in which the construction of a route is unsuitable. If these areas are forbidden then set N = +∞ . If these areas are not prohibited but unsuitable for construction or crossing, then N < ∞ is a sufficiently large positive number. For the computation of a shortest route regarding the construction of a start node v and destination node t, the following objective function has to be minimized:

3
Computing the trajectories for the development of optimal… for I i,j = 1, … , u i,j such that for all v ∈ V:

Cheapest routes
The construction of a route can depend on several factors, such as the construction price, the ground conditions on which the route is to be built, and personnel costs. In this paper, the weights of a grid ℝ 3 (U) = V, E, W C are the average construction costs W C of the edges E U . The cost of building a route per length unit is made up of the average construction cost per length unit and the ground condition. The objective is to compute a route between two nodes on ℝ 3 (U) that can be built most cost-effectively. The soil condition can be described by a set of discrete polygons P = P 1 , … P n on U where each P i = V i , E i consists of nodes K i and edges E i . Each polygon is assigned a color that describes the structure of the ground. For example, if the color of a polygon is j then this polygon represents a ground structure s j . All polygons of P with a color j create a polygon class G j ⊂ P . Each polygon of P is assigned to exactly one polygon class. All polygons of P are classified in k ≤ n different subsets so that P = G 1 ∪ … ∪ G k consists of k disjunctive polygon classes. Each polygon class G j is assigned a surface structure value j > 0 . If a route passes through G j , this value is the multiple of the construction costs of this route. This means if a route is built over an area P u of G j the construction costs on this area increase by the factor j . Let e i,l be an edge of s segments such that the edge e i,l = e i,l (1), … , e i,l s i,l can be decomposed into the segments e i,l (1), … , e i,l s i,l . The cost c i,l of an edge e i,l can be expressed as the sum of all incurred costs of the segments where K i,l = 1, … , s i,l is the index set of the segments, J i,l is the index set of the occurring surface structures on an edge e i,l , C is the average costs for the construction and maintenance of a route per length unit, d k is the length of the k-th segment of an edge e i,l . To identify a segment e i,l (k) on G j , an auxiliary node m k is generated in the middle of this segment. This point is checked by the Point-in-Polygon algorithm (Sunday 2001) to see if it lies in a closed polygon P u of P. If not, this procedure is continued for the next polygons until the result is positive. If m k is in P u with the color j, the segment e i,l (k) is assigned the value j . Calculating the length d k of a segment e i,l (k) lead to one of the following cases: 1. If there are no intersections and the auxiliary node m k is not in P u , then 2. If there are no intersections and the auxiliary node m k lies in P u , then the following is applied with k = 1 3. If there is exactly one intersection s 1 with e i,l (k) = v i , s 1 , then 4. If there is exactly one intersection s 1 with e i,l (k) = s 1 , v l , then 5. If there are two intersections s 1 and s 2 with e i,l (k) = s 1 , s 2 , then If a segment e i,l (k) lies in a polygon P u with the color j, then the segment length d k is multiplied by the corresponding surface structure value j . Note that the distances d k have to be calculated analogously to the equation (12). Figure 5 shows an edge e i,l lying on the three polygons P 1 in red, P 2 in green and P 3 in blue. It has two intersections s 1 and s 2 . The first intersection s 1 is at the boundary between P 1 and P 2 and the second intersection s 2 is between P 2 and P 3 . The edge e i,l is divided into three segments e i,l (1) , e i,l (2) and e i,l (3) . In the middle of these segments are the midpoints m 1 , m 2 and m 3 which determine the surface structure values 1 , 2 , and 3

Most convenient routes
In this part two variants for the calculation of the trajectory of a most convenient or passable route on U are presented. In the first variant a method is presented in which a route between two nodes is computed with a minimum sum of absolute elevation changes. This means that a route with the least vertical variations is searched for. Similarly, in the second variant, a route with a minimum sum of absolute gradients from the starting point to the destination is calculated.

Absolute elevation change
Both for ordinary traffic and for the transport of heavily loaded vehicles, routes are preferred which are easy and comfortable to drive on. Due to the topography of U, the surface of U has both ups and downs corresponding to the grid ℝ 3 (U) . An edge e i,j = e i,j (1), … , e i,j u i,j of V is divided into u i,j segments. Let V aux be a set of auxiliary nodes on the grid ℝ 3 (U) and A ⊂ (x, y, 0) ∈ ℝ 3 an area with an elevation of 0. Analogue to Sect. 3.1, the weight of a segment e i,j (k) is defined as an absolute elevation change between two neighboring nodes v i,j (k) and v i,j (k − 1) of V ∪ V aux as follows for k = 1, … , u i,j . A route between the starting node s and the destination node t is called the most convenient route regarding the absolute elevation change, if the sum of the absolute elevation changes on ℝ 3 (U) is minimal. Figure 6 shows that there are two nodes v i and v j on U and between these nodes there are horizontalequidistant auxiliary nodes v i,j (k) of V aux for k = 1, … , u i,j − 1 . These nodes generate between v i and v j the segments e i,j (1), … , e i,j u i,j on the surface of U. The surface of U is depicted between v i and v j by a red connection curve. An edge line consisting of the segments e i,j (1), … , e i,j u i,j approximates this connecting curve in order to discretely model the surface of U. Discrete modeling of U is sufficient since the 1 3 elevation data of terrain are also discrete. The red arrows illustrate the elevation changes between all neighboring auxiliary nodes along an edge e i,j .
If an absolute elevation change h i,j (k) of a segment of a route exceeds a maximum elevation change h max > 0 then this route may not be passable. A segment e i,j (k) is called critical with respect to the absolute elevation change if h i,j (k) > h max . An alternative most convenient route with h i,j (k) ≤ h max for all k ∈ 1, … , u i,j and e i,j ∈ E is searched where u i,j is the maximum number of segments of e i,j . All weights of the segments within ℝ 3 (U) will be modified when considering a maximum elevation change. Let h max > 0 be a maximal allowed absolute elevation change of a segment e i,j (k) . Then the weight of the segment e i,j (k) is defined as where N 1 , N 2 ∈ ℝ + are two sufficiently large positive numbers. The weight of an edge e i,j is given by If we assume that the absolute altitude change of e i,j (k) may not be greater than h max then N 1 = +∞ . If for e i,j (k) the areas from O are prohibited then N 2 = +∞ when this segment is in O or intersects this set. A route is searched between s and t of V with minimal absolute elevation changes: ∑ e i,j ∈E x e i,j w i,j ⟶ min! Fig. 6 A depiction of the decomposition of an edge e i,j = v i , v j in smaller segments with weights of absolute elevation changes. The red curve describes the topology U of a geographical area A, the black nodes v i,j (1), … , v i,j (u − 1) are the auxiliary nodes between v i and v j on U and the blue connections are the segments of an edge e i,j optimal solution defined. A set of solutions is determined in which the improvement of one objective function can be achieved only by the worsening of another, i.e., the set of optimal trade-offs (Pareto front). In our case, the Pareto front gives the optimum of the combination of the three criteria, the shortest route with respect to the distance, the cheapest route regarding the construction costs and the most convenient route in terms of the absolute elevation change or absolute gradient. In our investigation, we have defined a Pareto objective function f Pareto (E) on ℝ 3 (U) that consists of a distance function f D (E) , a cost function f C (E) , and an absolute elevation change function f H (E) as follows: with and such that for all v ∈ V: where + + = 1 for all , , ∈ [0, 1] . Since ℝ 3 (U) is discrete, the minimum of the objective function f Pareto (E) generates a solution space S Pareto of points in ℝ 3 for all , , ∈ [0, 1] with respect to the equation + + = 1 . This space lies on a hypersurface of ℝ 3 . Each point of S Pareto is a minimum of the function f Pareto (E) . The Cartesian coordinate axes that depicted the solution set S Pareto represent the three criteria distance, construction costs, and absolute elevation changes of a route. The values of each axes indicate the quantity or unit of the objective functions f D in meters or kilometers, f C in monetary units and f H in meters, respectively. The values of S Pareto indicate which values the respective objective functions have for each axis in ℝ 3 . In case one is looking for Pareto optima with two criteria with respect to {( , , ) ∈ [0, 1] 3 | + + = 1} then one can orthogonally project the set S Pareto from ℝ 3 to ℝ 2 . The projection of S Pareto is a curve in ℝ 2 in that each coordinate axis of this subspace indicates the values of the two remaining criteria. Let = c 0 ∈ [0, 1) be a constant with respect to the equation (23). Then = 1 − c 0 − for all ∈ 0, 1 − c 0 . The set S Pareto is projected for a fixed c 0 onto a curve in ℝ 2 . This will give all optimal compromises between the two functions f D and f C for any c 0 part of the function f H .

Computing optimal routes using CONTRA
For the construction of new routes, we developed a software called CONTRA , which is the acronym for Computing an Optimal Network of Transit Routes through mathematical Algorithms. This tool is capable of computing and planning the construction of new routes connecting two locations (e.g., cities) in any country. These routes consist possibly of different infrastructure segments such as roads, railways, gas pipelines, etc. This software uses GIS-data, e.g., ground conditions and elevation data, and optimization techniques to compute the building of optimal routes. The output of the software is a summary of the data of the computed route, such as road coordinates, construction costs, and lengths. As an example, we compute routes between the city of Khas Uruzgan and the city of Kabul . In this part we computed routes that are below the elevation of 3500 m. In practice, it is unusual to build routes in elevated areas, because the oxygen decreases with increasing altitude.

Input data
Contra uses for the computation different GIS-data of a chosen area. To determine an optimal route, the following input data is required by this program: 1. Coordinates of the start and destination points (in latitude and longitude) are given in Table 1. These two points are colored in red in Fig. 7. The red dashed line is the line-of-sight connection between them. 2. For the national land cover of Afghanistan the shapefiles of the organization AIMS (1997), formerly of the Afghan Geodesy and Cartography Head Office, were used (AIMS 1997). The whole country is partitioned in polygonal shapes that describe the respective type of land, see Fig. 8. 3. For the topographic representation of Afghanistan, the SRTM-30m data of NASA are used Watkins (2019). The resolution of the SRTM-30m data for Afghanistan is 1-arcsecond (in total 36001 × 54001 pixels) in a latitude/longitude projection (EPSG:4326). As an example, in Fig. 9 the area of the province of Uruzgan in Afghanistan is depicted. Computing the trajectories for the development of optimal… 2010). For our research we assumed an estimated cost depending on the different land covers. In general, it is possible to modify the estimated cost. 6. We assume the construction costs of roads along a mountain pass or on the mountainside on average USD 3000 per meter. The costs of quarrying and earthworks is higher than that of building an ordinary road.

Results
We compute optimal solutions for the construction of a route from the city of Khas Uruzgan to the city of Kabul in Afghanistan regarding three different single objectives: the shortest route in red, the cost-minimal or cheapest route in green, and the most convenient route w.r.t. the elevation change in blue, see Fig. 10. The columns in Table 2 describe (a) the construction costs of the routes in million USD, (b) the lengths of the routes in km and (c) the absolute elevation changes of the routes, i.e., the sum of all height changes from the starting point to the ending point along the route. Each of these three routes has a certain length and elevation profile. Figure 11 shows the elevation profile of a shortest road with respect to the distance between the cities Khas Uruzgan (left) and Kabul (right). Note that there is Fig. 9 The topography of the Uruzgan province in Afghanistan. The color spectrum from dark to light indicates the height of the terrain Fig. 10 Three single-objective optimal roads between the cities Khas Uruzgan and Kabul in Afghanistan computed with CONTRA : shortest (red), cheapest (green), and most convenient (blue). Image rendering: Google Earth  Table 2) and the re-rendering of that very path in Google Earth (299 km). We believe that this is due to the fact that Google Earth does not take the elevation changes into account when computing a path length. Figure 12 shows the elevation profile of a cheapest road with respect to the construction costs from Khas Uruzgan to Kabul. Analogously, Fig. 13 visualizes the elevation profile of the most convenient road above sea level between these two cities. One can see that this route is a longer distance, but most parts of the route have low absolute elevation changes. For the Pareto optimum, the following objective function of cost and elevation change is minimized for different 0 ≤ ≤ 1: Figure 14 shows the Pareto front of the routes from Khas Uruzgan to Kabul regarding the two objectives "construction costs" and "elevation changes". This chart shows on the horizontal axis the construction costs of the routes in million USD and on the vertical axis the absolute elevation changes in meters. The small red circles represent specific routes between these two cities. The leftmost circle represents the cheapest route and the right most circle represents the most convenient route with minimal altitude variation. The circles that lie between these two extremal circles, are other optimal routes from Khas Uruzgan to Kabul which are trade-offs of these two criteria. The gaps between the circles in this figure indicate that there is an unsuitable area (obstacle) between two routes so that the construction of a route across this area is sub-optimal. Thus, the routes run around such "obstacle areas" and generate gaps in the Pareto front. x e i,j c i,j + (1 − )w i,j ⟶ min! Fig. 12 The elevation of a cheapest road on terrain from left Khas Uruzgan to right Kabul in Afghanistan. Image rendering: Google Earth Fig. 13 The elevation of a most convenient road with minimal absolute elevation changes from left Khas Uruzgan to right Kabul. Image rendering: Google Earth Figure 15 shows the optimal routes picked from the red circles in Fig. 14 for six different values. For the extremal value = 1 the Pareto optimal route is the cheapest route between Khas Uruzgan and Kabul (see top left depiction in Fig. 15). The second extremal value = 0 presents the most convenient route between these two cities (see below right depiction in Fig. 15). All the other depictions are optimal compromises between these two extremes.

Conclusions
The search for a suitable trajectory for the construction of a route in the transport or energy sector can be very time-consuming considering the environmental information. This paper presented a new method of computing the shape of a new route using discrete optimization and mathematical algorithms. The relevant information was taken from the GIS data such as the geographical environment, structure of the soil, topographic elevation data of terrain, and average construction costs of an infrastructure segment. In the first step, U was discretely approximated by a weighted grid . In the second step, the weights of this grid were computed regarding the three criteria: (a) distance, (b) construction costs and (c) absolute elevation changes or absolute gradients. In the third step, routes with minimal edge weights were determined between two locations whereby these routes are globally optimal. Furthermore, with Pareto optimization, we have Absolute elevation change in m

Pareto curve for optimal routes
Fig. 14 The Pareto front of the most convenient and cheapest routes regarding absolute elevation change and construction costs from Khas Uruzgan to Kabul. Each red circle represents an optimal trade-off routes between the cheapest and the most convenient routes where the most left circle represents the cheapest route and the most right circle represents the most convenient route combined these optimization criteria to determine optimal compromises between them that correspond to the real conditions. Finally, we modeled the area between Khas Uruzgan and Kabul in Afghanistan with our computer program CONTRA and with the corresponding GIS data and determined the shortest, cheapest, and most convenient roads between these two cities. With the Pareto optimization we determined several compromises between the cheapest and the most convenient route and presented the change between both extremes as trajectories. In further work, we have the goal to compute the construction of optimal routes with tunnels and bridges. The calculation of tunnels or bridges along a route is more complicated than routes on the surface. When planning tunnels and bridges, the