Advancing Transportation Routing Decisions Using Riemannian Manifold Surfaces

We consider several real-world driving factors such as the time spent at traﬃc signs (e.g., yield signs and stop signs), speed limits, and the topology of the surface to develop realistic and accurate routing solutions. Though these factors increase the complexity of modeling, they provide the ﬂexibility to evaluate the routing solutions from diﬀerent perspectives: cost, distance, and time, to name a few. First, we develop a set of algorithms based on the Riemannian manifold surface (RMS) to factor in the Earth’s curvature to calculate distances. Second, we present a multiobjective, nonlinear, mixed-integer model (MINLP) that minimizes the distance traveled, time traveled, traveling costs, and time spent on traﬃc signs to design and evaluate the routes where the waiting times associated with traﬃc lights, stop signs, and yield signs are stochastic. Finally, we apply MINLP and RMS-based algorithms to a set of real-life and short- and long-distance transportation problems and analyze the results from computational experiments and discrete event simulations. We show that our approaches are on par with the state-of-the-art application, Google Maps, and yield realistic routing solutions that generate signiﬁcant cost savings.


Introduction
Transportation is one of the essential activities in a supply chain. For example, goods manufactured at plant locations must be distributed to customer locations through warehouses, distribution centers, or cross-docking facilities. Typically, a fleet of trucks (owned or outsourced) performs the task of distributing goods. In the view of fuel prices, the complexity of the supply chain, and the uncertainties in the market place, the trucking operations must be efficient in order to keep costs under control. For instance, in 2019, US companies spent 668 billion dollars on trucking costs alone [1]. As a result, businesses are continually looking for a new approach or technique to reduce transportation costs. In product distribution, the sequence (or route) in which goods are delivered by a truck determines the distance traveled by the truck, which, in turn, determines the cost of transportation. e distance between locations is often used as a surrogate measure for the transportation cost. It is therefore vital to understand how distances are calculated. e assumptions related to distance calculations play a key role in determining the accuracy of the measurement. One of the assumptions is the surface on which the distance is measured. Often, the Euclidean surface is used because it is less complex to calculate distances [2]. However, the Euclidean surface generates inaccuracies in distance calculations because it does not take into account the curvature of the Earth. As a result, the use of the Euclidean surface has implications such as suboptimal routing decisions, resulting in increased transportation costs. Motivated by the above, we develop an approach based on the Riemannian manifold surface (RMS) that factors in the curvature of the Earth when calculating distances. For instance, the curvature of a Riemannian manifold surface indicates how much Earth's surface is curved in a local neighborhood. As the surface changes, the curvature changes. We use this curvature change to calculate the shortest path distances.
Subsequently, we demonstrate that the RMS-based approach improves the distance calculation, yields realistic distances, and generates significant transportation cost savings.
In this work, to factor in the curvature of the Earth, we use the generalized notion of a "straight line" to "curved spaces," i.e., a geodesic. e calculation of geodesic distances involves the use of local topological regions on Earth's surface. For example, while sprinting on multilane running tracks, a runner on the inner lane starts at the back, while a runner on the outer lane starts at the front. e starting points of the runners are determined by factoring in the geodesic distances, which are two dimensional as the running tracks are flat. However, driving a vehicle involves three-dimensional geodesic calculations that we use in this research.
Although companies use distance as a critical measure to evaluate the routing decisions, end consumers use various indicators (criteria) to determine routes for work and leisure travel. For instance, travel time, travel costs, travel distance, and waiting time spent at traffic signs (e.g., stop signs, yield signs, and traffic lights) are often used. Nevertheless, existing global positioning and global information systems do not provide, for example, the impact of traffic signals on travel cost and time and the flexibility to evaluate routing decisions using various criteria. Moreover, these systems do not consider local shortest geodesic distances. Inspired by the above, we develop a multiobjective, nonlinear, mixed-integer model (MINLP) where the end users can evaluate the routes with multiple criteria based on travel distance, cost, time, and traffic controls such as stop signs, traffic lights, and yield signs. Next, we illustrate the RMS-based approach and the MINLP model, and later, we perform simulation experiments using real-life instances.

An Illustration of Riemannian Manifold Surfaces and
Multiple Criteria Route Evaluation. In this section, we illustrate two key concepts that are used in our methodology: (i) Riemannian manifold surface-(RMS-) based approach to calculate the distance; (ii) multiple criteria to evaluate the routes. For illustration purposes, we choose two points, Origin, O, and Destination, D. Let O and D be Hamden and New Haven, two cities in the state of Connecticut, for the ease of illustration as the authors are familiar with the geography. ere are three routes between Hamden and New Haven: (A) via I-91 South-using highway; (B) via Whitney Avenue-using local roads; (C) via Ridge Road-using local roads.
To illustrate the concept of RMS, we focus on Route A. Figure 1 is a schematic representation of Route A. In this figure, the solid lines represent the boundaries of the road, and the road is divided into two lanes: Lane 1 and Lane 2. For instance, a vehicle traveling from the origin to the destination can either start from Lane 1 or 2. If the vehicle starts from Lane 1, switches to Lane 2 when the road curves, and stays in Lane 2 until it reaches the destination, the path traced by the vehicle is referred as the innermost lane. Conversely, if the vehicle starts at Lane 2, switches to Lane 1 when the road curves, and stays in Lane 1 until it reaches the destination, the path traced by the vehicle is referred as the outermost lane. e dashed lines in Figure 1 represent the inner-and outermost lanes of the road. e distance traveled by the vehicle using inner-and outermost lanes is 10.84 miles and 11.38 miles, respectively. Note that both the inner and outer lanes are created based on the shortest and the longest geodesic paths, respectively, that takes the Earth's curvature into account which is a central theme of the RMSbased approach (the details are provided later in the paper). e approach based on RMS yields savings of 11.38 − 10.84 � 0.53 miles for Route A. Similarly, savings on Routes B and C are identified using RMS, and Table 1 summarizes the results. Note that when two locations are only a few miles apart, the impact of the Earth's curvature on the distance calculated is very minimal. is can be seen in 1, as the difference between the outermost and innermost lanes is less than a mile for all routes. Later, in the numerical section, we show that the difference between the outermost and innermost lanes (i.e., distance savings) can be significant for long-distance routes. In such cases, it is feasible for a traveler to devise a driving plan to realize the savings in the distance.
Next, using the above real-world instance, we illustrate how a decision maker can use multiple criteria to evaluate routes A, B, and C. We consider the following criteria: distance traveled, waiting time at traffic signs, travel cost, the number of traffic signs, and travel time. Table 2 provides a summary of the evaluation of the routes, and the associated details are in Section 6. Note that for each criterion, the best values are marked in bold.
For instance, if the distance traveled or the traveling cost is used as a sole criterion, then Route B is the best. However, Route B is the worst when the time traveled is considered. In this case, both A and C are better than B. When the waiting time at traffic signs is considered, Route C is the best. For instance, Route A has 1 stop sign, 1 yield sign, and 16 traffic lights which totals to 18 traffic signs. Route B and C has a total of 29 and 21 traffic signs, respectively. If the number of traffic signs is used as a criterion, Route A is the best. e above illustration highlights the importance of using multiple criteria in route evaluation as the route selection varies with the criteria considerably.  In summary, we develop a set of algorithms that are based on Riemannian manifold surfaces to address the need to calculate the distance accurately. In addition, we formulate a mathematical programming model that considers the real-world driving factors such as the time spent at traffic signs (e.g., yield signs and stop signs) and speed limits to evaluate routing decisions from multiple perspectives: cost, distance, and time, to name a few. Furthermore, we apply our approaches to real-life transportation problems and analyze the results from computational experiments. e rest of this paper is organized as follows: Section 2 summarizes the relevant literature and the research contributions of this paper. Sections 3 and 4 offer details on the RMS-based approach to calculate the distance. Section 5 provides details on mathematical programming formulations that use multiple criteria to evaluate the routes. Section 6 summarizes computational experiments and simulation results. Section 7 concludes with a general discussion of managerial insights and avenues for future research.

Literature Review
In this section, we first present the literature related to the distance calculation since one of the main contributions of our research is in developing a new methodology to calculate distances. Next, we look into relevant literature where multiple criteria are used to evaluate routing solutions. Finally, we present our contributions in relation to the past work.

Distance Calculation Methodologies.
In practice, there have been several software applications developed to calculate distances using travel information such as oncoming traffic situation, resting areas, and time to destination. However, the visibility into mathematical methodologies that are used in such applications is very minimal. erefore, it is a challenging task to subject the methodologies used by software to rigorous academic discussion. On the contrary, in the academia, research related to distance calculation methodologies is limited. In the papers where the distance calculation is addressed, typically, rectangular or Euclidean distance functions are used. In their seminal work, Love and Morris [3] provide mathematical proofs for calculating transportation distances using convex functions that extend beyond the rectangular or Euclidean distance calculations. Brimberg and Love [4] consider a new distance measure, the weighted one-two norm, which is a positive linear combination of the preceding norms as compared to the Euclidean or rectangular distance calculations. Fernández et al. [5] compare the weighted l p -norm, which has been proven to be an accurate distance calculation function, to the l 2b -norm which uses different parameter estimates for the distance prediction. is methodology is extended by Brimberg and Love [4] andÜster and Love [6]. Other papers such as Von Hohenbalken and West [7]; Miyagawa [8]; Alpaydin et al. [9]; Griffith et al. [10]; Herica and Servio [11]; Havelock et al. [12] have all used different methodologies to calculate distances, from the Manhattan and Euclidean distances to Lobachevskian and p-norm methodologies.

Multiple Objective Routing
Solutions. An extensive amount of research has been done in formulating methodologies to aid the decision maker with optimal routing decisions; see, Toth and Vigo [13] for a review on vehicle routing problems (VRP). We review only the research that considers multiple criteria/objectives to evaluate the routing solutions. For instance, Bowerman et al. [14] develop an approach to model the urban school bus routing problem that considers multiple criteria such as route length, number of routes, load balancing, length balancing, and student walking distance. Zografos and Androutsopoulos [15] identify routing solutions that minimize both transportation time and risk as a part of a decision support system to transport hazardous materials and address emergency situations. Demir et al. [16] consider a variant of the pollution routing problem that minimizes both the fuel consumption and driving time. Kovacs et al. [17] extend the VRP from a service perspective by considering several objectives related to improving driver consistency, arrival time consistency, and cost minimization. Kumar et al. [18] consider a production and pollution routing problem with time windows that minimizes the total operational cost and emissions. See Jozefowiez et al. [19] for an extensive review on research related to multiobjective optimization in routing problems. ough our research uses criteria that are related to the total cost and the distance, it also considers the number of traffic signs and the time spent on these signs as criteria that have not been considered in the past research.  As mentioned above, one of our primary goals is to highlight the lack of research in using criteria related to traffic signs and to stress the significance of using multiple criteria when making routing decisions. erefore, we restrict our attention to modeling multiple criteria and demonstrating the model's usefulness by applying our techniques to a few real-life instances. We note that our research does not involve designing heuristics for the routing problem. As future research, our work can be extended by employing approaches such as genetic algorithmbased techniques (e.g., [20,21]), greedy heuristics with variable neighborhood search (e.g., [22]), learning-based heuristics (e.g., [23]), and fuzzy programming (e.g., [24]), to name a few, to solve large-scale, multiobjective, routing problems.

Research Contributions of
is Paper. Our work is closely related to the study by Alpaydin et al. [9], Griffith et al. [10], and Havelock et al. [12]. Alpaydin et al. [9] propose a nonparametric approach using neural networks for estimating actual distances compared to the known distance estimators which are parameterized functions of the coordinates of the points. Griffith et al. [10] explore the non-Euclidean nature of network spaces, including Lobachevskian space (via the Poncaŕe disc), and compare the Euclidean, Manhattan, Minkowskian, and Lobachevskian distance calculation predicating parameters. e authors conclude that the skewed transportation networks in largescale regional landscapes are characterized by a Lobachevskian geometry, whereas most smaller scale landscapes appear to be best characterized by a Minkowskian space whose parameters are between those of a Manhattan and Euclidean space. Havelock et al. [12] propose the use of learning-based strategies to solve the problem of distance calculation by using a stochastic method called the adaptive tertiary search (ATS) strategy. is methodology works by utilizing the information provided in the coordinates of the nodes to calculate the actual distances.
While we build on the previous work, our research differs in the following ways: (1) We develop a novel methodology based on the Riemannian manifold surface (RMS), a concept that uses the road network topology and factors Earth's curvature into account. As part of our methodology, we propose a set of algorithms where we adapt and modify the linked chain methodology (LCM) developed by Tokgöz et al. [25] to calculate distances. (2) We formulate a multiobjective, nonlinear, mixedinteger model (MINLP) to design and evaluate the routes where distance traveled, time traveled, traveling costs, and time spent on traffic signs are minimized. In this model, the waiting times associated with traffic lights, stop signs, and yield signs are stochastic. (3) We apply our approaches to a set of real-life and short-and long-distance transportation problems and analyze the results from computational experiments and discrete event simulations. We show that our approach yields realistic solutions and generates significant cost savings compared to using the state-of-the-art applications (e.g., Google Maps).

Riemannian Manifold Surfaces and Linked Chain Method
Carl Friedrich Gauss proved eorema Egregium (meaning remarkable theorem in Latin) in 1828 by establishing an important property of surfaces. e theorem states that the Gaussian curvature of a surface can be determined entirely by measuring distances along the paths on the surface, and the curvature does not depend on how the surface is embedded in a three-dimensional space. Later, Bernhard Riemann extended Gaussian theory to higher-dimensional spaces that allows distances and angles to be measured and the notion of curvature to be defined. e Riemannian approach is an intrinsic way to the manifold (a topological space that is locally Euclidean) and does not depend upon how the manifold is embedded in higher-dimensional spaces. A Riemann manifold M is a metric space on a smooth manifold with a smooth section of the positivedefinite quadratic forms on the tangent bundle when the length of a continuously differentiable curve c: [a, b] ⟶ M is defined by the following equation: e definition of length in (1) allows every connected Riemann manifold M to become a metric space on which the distance between a and b on M is defined as the infimum of l(c) when c is a continuously differentiable curve joining a and b. Later, we use (1) to calculate the geodesic length which is instrumental in calculating route distances. One of the most popular use of the Riemannian manifold is in the development of general theory of relativity by Albert Einstein; in particular, his equations for gravitation are constraints on the curvature of space-time.
e curvature of a Riemannian manifold surface (RMS) indicates how much the surface is curved in a local neighborhood. For instance, the Riemann manifold is a plane in the Euclidean space when the curvature is zero; however, when the curvature changes from zero to one, the surface changes from a plane to a 3-dimensional sphere. is change in the curvature is important when calculating shortest path distances on Earth's surface. is surface is a natural manifold, and the pathways on this manifold are natural geodesics. e research on RMS-based approaches in transportation is almost nonexistent except by Tokgöz et al. [25]. For instance, Tokgöz et al. [25] solved a location routing problem (LRP) on RMS and showed that their solution methodology generated results that are more realistic than using Euclidean surfaces since the Earth's surface has changing local RMS curvatures. e authors developed an approach called the linked chain method (LCM) to calculate geodesic distances. We adapt and modify the LCM approach to our problem. Typically, in the LCM approach, given an origin (O) and destination (D), circles are generated to cover the focus area that contains the O-D pair. Unlike Tokgöz et al. [25], these circles are not randomly generated, and in addition, the circles represent the underlying road structure in our case. For instance, using the LCM, we generate innerand outermost lanes for highways that provide a range of traveling distance which has not been addressed in the past research. Subsequently, the circles that cover geodesics are identified, and the length of the geodesics is calculated. e overlapping circles that contain geodesics are identified, and they are linked (hence, this approach is called the "linked chain") to create a path that connects O and D. Next, we provide details on the linked chain method.

Algorithmic Solution for Shortest Path Geodesic Distance Calculation
In this section, we introduce four algorithms to determine the shortest path between two points: origin, O, and destination, D, that can be used for routing vehicles. In addition, there are two theorems and a corollary stated and proven by using isomorphism and LCM-related routing results. e set of notations used in the development of algorithms is given in Table 3. We assume geodesics are determined based on the pathway intersections and illustrate this assumption by using a T-intersection. For instance, a T-intersection, see Figure 2, has three arms, and each arm could be a geodesic. Note that geodesics can be curved lines, but for illustration purposes, geodesics are shown as straight lines, see Figure 2. We assume that there exist nonoverlapping circles containing all the geodesics; these circles form the topology on the manifold. As a result, each circle can have only one geodesic. For instance, in Figure 2, three geodesics (three arms) are covered by three different circles. e radius of these circles is different due to varying lengths of the roads that form the T-intersection. Moreover, if a path uses the T-intersection, two circles of varying radii that cover two geodesics will be a part of this path. For instance, geodesics 1 and 2 could be one path and geodesics 1 and 3 can be another. Our first objective is to minimize the total distance traveled which is calculated using (1) by replacing a and b with O and D; i.e., 4.1. Subroutine 1. Algorithm 1 shows the pseudocode for the first subroutine. First, a surrounding area M 0 within M of O and D is selected. Next, r is set to 0.1, and circles with radius r are drawn so that the entire surrounding area is covered with circles. Some of these circles may have geodesics in them and some may not. Once geodesics are identified for a radius r, i.e., c ir 's, they are added to a set G. is process of identifying geodesics is repeated by incrementing r, i.e., r � r + ϵ where ϵ is set to 0.001. Note that a geodesic that is identified by multiple radii will be added only once to G. As a result, geodesics in G are not duplicated. Each one of these geodesics could have a minimum and a maximum length based on its curvature, and these lengths are calculated using Algorithm 4. Figure 3 illustrates the steps involved in Algorithm 1. For illustration purposes, a focus area M 0 that is similar to a busy neighborhood in Manhattan, New York, is selected. In this illustration, the focus area is represented by a rectangle that has 8 intersections: four on the top row and four on the bottom row, see Figure 3(a). First, the entire M 0 is covered with the circles of initial radius r 0 . For the purposes of   Figure 3(a) is not shown; therefore, some of the areas in the figure are not covered. Next, the circles with geodesics are identified. In Figure 3(a), four geodesics are identified (on the top row) by four circles that are indicated by the dashed lines. We assign circles C c 1r 0 , C c 2r 0 , C c 3r 0 , and C c 4r 0 to geodesics c 1r 0 , c 2r 0 , c 3r 0 , and c 4r 0 , respectively. Finally, the above geodesics are added to set G. e above steps are repeated with a new radius, r 0 � r 0 + ϵ. If geodesics are identified with the circles of new radii, they will be added to G, otherwise, not. Figure 3(b) illustrates the geodesics that are identified when the radius is equal to r 1 . In this figure, the circles that identify geodesics are represented by the dashed lines. We assign circles C c 1r 1 , C c 2r 1 , C c 3r 1 , and C c 4r 1 to geodesics c 1r 1 , c 2r 1 , c 3r 1 , and c 4r 1 , respectively. We add the above geodesics to set G.

Subroutine 2.
Algorithm 2 shows the pseudocode for the second subroutine. We use this subroutine to calculate the length of geodesics, c, that are in G. Initially, two empty sets are created, G new and P. First, in set G, we select a radius r and examine the geodesics, c ir , identified with that radius. Second, the length, l(c ir ), of every geodesic is calculated using the following function: ird, we define p c ir to be a vector with three elements: the geodesic c ir , the length of the geodesic l(c ir ), and the circle C c ir that covers the geodesic. Subsequently, we add the geodesic c ir to set G new and delete the same geodesic from set G. Finally, p c ir is added to set P. e steps explained above are repeated until the distances for all the geodesics are calculated, i.e., until G is empty. As an output from this subroutine, we identify a set P that consists of vectors that contain the geodesics, length of geodesics, and unique circles that contain the geodesics.

Subroutine 3.
In topology, geodesics are connected and form a path if circles that contain them are overlapping.
is subroutine searches for overlapping circles to determine the paths within the focus area, M 0 . Algorithm 3 shows the pseudocode for the subroutine. Initially, an empty set, N, is created to store all possible paths, i.e., linked list vectors of geodesics. en, two elements p c ir and p c jr of set P are chosen such that i ≠ j. If circles C ir and C jr that contain c ir and c jr are overlapping, the link is added to set N; otherwise, a different pair of p is chosen. is process is repeated until all pairwise comparisons in P are exhausted. As a final output, a set N that consists of all possible paths is generated. (1) for r � r 0 to max|M 0 | do (2) Cover M 0 with LCM-linked circles with radius r by tracking all possible roads between O and D (3) for i � 1 to n r do (4) Assign c ir to be the i th geodesic covered for radius r (5) Assign C c ir to be the circle that covers c ir (6) if c ir ∉ G then end if (11) end for (12) end for ALGORITHM 1: Geodesic search pseudocode. 6 Journal of Advanced Transportation all c ir from intersection to intersection which are then used to calculate the minimum and maximum length route between O and D. Algorithms 1-4 yield the following set of theorems.
Theorem 1. C and G are isomorphic.
Proof. ere exists a set of unique circles, C, that covers the set of geodesics within the given bounded region. One circle can cover only one geodesic, and one geodesic can be covered by only one circle. It is evident from Algorithm 2 that there exists a unique circle per geodesic; therefore, there exists a map f : C ⟶ G that is one-to-one. To show the converse, assume there are two circles in G that correspond to the same geodesic. In this case, two circles must have the same overlapping intersections which would conflict with Algorithm 1 implementation. erefore, f is an onto map; there exists a 1-1 and onto map that projects C to G. e (1) for r � r 0 to max|M 0 | do (2) for i � 1 to n r do (3) for c ir in G do (4) Calculate l(c ir ) end for (10) end for (11)  □ Theorem 2. LCM determines the distance minimizing shortest path geodesic between two points chosen on a compact path-connected manifold region on Earth's surface.

Proof.
A compact path-connected manifold region on Earth's surface indicates the existence of a path that allows the existence of geodesics between the two points chosen in a compact region. Algorithm 2 is used to determine all possible pathways between the two points, while Algorithm 4 is used to determine the distance minimizing shortest path on the compact pathconnected predetermined region. e set N h in Algorithm 4 is the set of all pathways on the manifold's surface.
□ Corollary 1. LCM applied on Earth's surface does not always yield the shortest path geodesics between the two points.
Proof. Assume there exists a geodesic that yields the shortest path between the two chosen points on Earth's surface. Such a geodesic would require to identify predetermined pathways with finite lengths and nonoverlapping loops; however, the looping nature of the manifold region could result in infinite distances. For instance, consider a region that only contains ocean. In this case, the chosen manifold's curvature is always 1; therefore, there exist infinitely many geodesics on this surface due to the looping nature on a spherical surface. Such a compact region of Earth's surface can generate infinite geodesics and infinite pathways. In this case, Algorithms 1-4 are not applicable.

A Multiobjective, Mixed-Integer, Nonlinear Mathematical Model
In this section, in addition to (2), we introduce four objective functions that minimizes (i) total time spent on traffic signs, (ii) total travel cost, (iii) total number of traffic signs, and (iv) total travel time. As the additional objective functions Define if D ∈ C c ir then (10) Repeat steps 3-6 by replacing O with D Find the minimum and maximum of l(c ir ) calculated in Algorithm 2 (15) Calculate the total length of the geodesics for all N h using step 14 (16) l min � min N h and l max � max N h are the minimum and maximum length geodesics for all N h that connect O and D (17) end if (18) end for ALGORITHM 4: Path finder and distance minimizer pseudocode.

Run: Algorithm 2
Input N � ∅ (1) for i � 1 to n r do (2) while p c ir ∈ P do (3) for j � 1 to n r do (4) if i ≠ j and C c ir ∩ C c jr ≠ ∅ then (5) Link c ir and c jr and add the link to N (6) Continue updating N with links for p c ir ∈ P (7) end if (8) end for (9) end while (10) end for ALGORITHM 3: Geodesic path construction pseudocode. 8 Journal of Advanced Transportation involve traffic signs such as stop signs, traffic lights, yield signs, and the speed limits, we introduce the assumptions and parameters that are associated with them. In some cases, objective functions such as travel distance and time might be correlated, but having them as individual functions provide the flexibility for the decision maker. For instance, in the densely populated areas (e.g., northeastern parts of US), mostly travel time and distance are not correlated. Note that we do not model road congestion in our work which could be an avenue for future research; however, the impact of the network congestion is captured by the objective functions such as travel time and cost. e set of notations used in the mathematical model is given in Table 4. e expected duration of a stop sign cannot be zero for lawful reasons; therefore, we assume that the wait time of a vehicle on the j-way stop sign, w j , a stochastic parameter, varies between waiting times min w j and max w j . e set of all stop signs is

Assumptions and Parameters
where represents the set of j − way stop signs k u j−1 on route u for all 1 ≤ u ≤ b and 1 ≤ j ≤ 4. For instance, a specific choice of a vehicle's wait time on a j-way stop sign can vary between waiting times min w j � 3 and max w j � 3j + 3 seconds. (ii) Traffic Lights. As one may have to wait at a traffic light longer than another (e.g., busy vs. nonbusy intersections), the maximum time spent can vary among traffic lights. Hence, it is reasonable to use the maximum time spent at traffic lights to categorize them. We categorize the traffic lights based on their maximum durations. Let the set represent the category-i traffic lights based on β u i and the maximum amount of time that a vehicle can spend waiting on red-light duration on route u such that 1 ≤ u ≤ b. For instance, Δ 100 1 represents the 1 st category of traffic lights that has a maximum duration of 25 seconds on the 100 th route. e set of traffic lights for b different paths from O to D takes the form for all 1 ≤ u ≤ b and e ∈ I 0 where I 0 is the number of traffic light categories. e number of traffic lights on the u th pathway in category i is represented by |Δ u i |, and the total number of traffic lights on each path u from O to D is for all 1 ≤ u ≤ b. For the example given above, if the set Δ 100 1 has 5 elements, then |Δ u 1 | � 5 indicates that there are 5 traffic lights with 25 seconds of maximum wait time on route 100. Let t u i represent the stochastic parameter that varies between 0 and β u i for waiting at the category-i traffic light on route u. (iii) Yield Signs.
ere are two types of traffic yields considered for a vehicle to join traffic: the yield signs to enter the highway traffic and the city traffic. Let y u 1 represent the number of yield signs to enter the highway and y u 2 represent the number of yield signs to enter the city traffic on route u for all 1 ≤ u ≤ b. Let y u � y u 1 + y u 2 for each u, and let μ u y k be the stochastic parameter representing the wait time at y u k for k � 1, 2. For instance, the average time it takes for a vehicle to enter a highway when there is a yield sign can be assumed to vary between 4 and 8 seconds, and the time it takes to enter the city traffic when there is a yield can be assumed to vary between 5 and 10 seconds on average. (iv) Speed Limit. e distance calculations depend on speed limit, geodesic distance, and traffic sign durations. e LCM (linked chain method) circles are used to determine the geodesic distances from one intersection to the other on each route. Each LCM circle contains the corresponding road's speed limit. We assume the set of speed limits to be where there are q different number of speed limits that could occur from O to D. If there are more than one in a LCM circle, then we consider an average of speed limits for time-related calculations. For instance, using geodesics c s from Algorithm 2 (subscript ir is modified as s) and the corresponding length l(c s ), the distance traveled on the LCM circle C c s can be calculated as follows: (l(c s ))/v l , where v l ∈ V is the speed limit value on the road covered by C c s . Later, we use this formula in the MIP formulations to determine the travel time from O to D. be the fixed per-mile fuel cost and per-mile operational cost, respectively. e constrained, mixed-integer, nonlinear program formulation for the b-path problem to travel from O to D that includes z 1 can be formulated as follows:
Recall that algorithms from Section 4 are used to determine the shortest geodesic path, and the corresponding objective function z 1 minimizes the travel distance. We reintroduce z 1 in (10) for clarity. In this formulation, u is an integer variable that denotes the route. In addition, c(t) in (10) is a nonlinear, continuous variable that is used to calculate the geodesic distance on nonlinear curves. e waiting times associated with traffic lights, stop signs, and yield signs are stochastic parameters. Since the number of traffic signs on each route is well-known and do not vary, they are constant values. Equation

Weighted Sum Formulation (WSF).
e weighted sum formulation (WSF) allows the decision maker to choose a route based on the weights assigned to the objective functions z 1 through z 5 . Function f(.) introduced below normalizes the objective function values z 1 -z 5 for each route, and it is used as a part of minimizing the weighted sum. is option of formulation requires weighted sum formulation of equations (15)- (19) by introducing the weight function σ h : A 1 ⟶ R + with the corresponding objective function: where A 1 is the finite index set and f(z u h ) are the normalized values of the objective functions z u h calculated for each route u by using the formula: for each h � 1, . . . , 5 and each u such that 1 ≤ u ≤ b. e values of f(z u h ) are determined for all h and u that forms a set; this set of values is used for determining the minimum value of z 6 . e weights σ (1) h are known constants that are assigned by the decision maker, and the sum of the weights must be equal to 100%.

Best Vector Set Solution (BV).
In this section, we propose a best vector set solution approach where the vector set consists of the best solutions determined by minimizing each z i , ∀i � 1, . . . , 5. Since the optimal route can vary with z i , the corresponding vector set solution can be different. In this case, where i, j � 1, . . . , 5, u � 1, . . . , b, b is the total number of the routes, and the best vector set BV � min(z u 7 ) which is a subset of all solutions. For instance, assume b � 4. Let route 1 (i.e., u � 1) be optimal for z 1 and z 2 , route 2 (u � 2) be optimal for z 3 and z 5 , route 3 (u � 3) be optimal for z 3 , z 4 , and z 5 , and route 4 (u � 4) be not optimal for any z i . e best vector (BV) solution is given as follows: Note in BV that there is no minimum objective function value for route 4 as it is not optimal for any z i . e routes u � 1, 2, 3 in the BV solution may or may not be displayed in the corresponding vectors listed in the set. Based on the end user's preference, a specific route could be picked from the vector set. For instance, if the user is interested in saving time and money, a route (s) that minimizes both travel time and travel cost can be chosen from the vector set. Later, we illustrate the BV solution for two real-life routing problems in Section 6.

An Algorithmic Solution.
In this section, we develop an algorithmic solution for the mathematical models presented in Section 5.2. Recall that equations (10), (12), and (14) are nonlinear since the distance calculations of geodesics depend on nonlinear curves. e numbers of traffic controls such as stop signs, yield signs, and traffic lights are known constants for each route. erefore, the objective function in equation (13) is a constant for a route. Algorithm 5 uses the outputs of Algorithms 1-4 from Section 4.
In this algorithm, constraints (15)-(18) introduced in Section 5.2 are randomized by using the minimum and maximum values (that are constant values and assumed to be known) for traffic signs. ese constraint values are then used for calculating z 2 −z 5 . Values of σ h are the inputs to the algorithm for calculating the weighted sum value z 6 . e values attained from z 1 −z 5 are then used for calculating the best vector set solution z 7 ; the decision maker can choose a route from a variety of routes by using z 7 that depends on the minimum values attained for z 1 −z 5 . A real-life application of Algorithm 5 with specific entries for minimum and maximum values is given in Section 6.
For instance, the set P which is an output of Algorithm 2 is a critical input to Algorithm 5. e set P consists of p c s vectors which has the following information: (i) the geodesic c s , (ii) the length of the geodesic l(c s ), and (iii) the LCM circle C c s that has c s . e stochastic parameters are calculated by using a uniform distribution.
Note that running Algorithm 5 for sufficient number of times and averaging the outcomes will yield the desired objective function values. e highest and lowest values of the parameters presented in constraints (15)-(18) determine the minimum and maximum values for the objective function z 2 . e locations and numbers of traffic signs along with the speed limits and maximum durations of traffic lights are assumed to be known.

Numerical Analysis and Simulation
In this section, we apply the methodologies from Sections 4 and 5 to two different real-life routing cases, short and long distances, and analyze the results from computational experiments and discrete event simulations. e first case includes three different short-distance routes that connect a location in Hamden (O) and another location in New Haven (D), two small cities located in the state of Connecticut. One of the routes is via a highway, while the other two routes are via local roads. For this O-D pair, we assume a travel pattern of two trips per day for all year. e second case includes three different long-distance routes that connect an origin and a destination in the states of Connecticut (CT) and California (CA), respectively. Note that the set of long-distance routes is chosen to mimic the actual truckload routes of a multinational stainless steel and specialty metal manufacturer located in the state of Connecticut. For this O-D pair, we assume a travel pattern of one-round trip per two weeks for all year. In both shortand long-distance routing cases, though there are many routes existing between O-D pairs, we restrict our attention to only three routes to demonstrate the theoretical results and the benchmark with the state-of-the-art application, Google Maps ™ . A summary of parameter values that are used in the computational experiments (unless otherwise stated, the parameter values are identical for both short-and long-distance routing cases) is given in Table 5.
Several types of traffic signals are present in each route: traffic lights (TL), yield signs (YS), and stop signs (SS). Furthermore, under the SS category, four different types of stop signs are considered. e time spent in waiting at the traffic signals varies with the type. We assume minimum and maximum waiting times for each type of the traffic signals: traffic lights (15 seconds, 75 seconds), yield signs (3 seconds, 9 seconds), and j-way stop signs (3j seconds, 3j + 3 seconds) where j is a type of a stop sign for all j � 1, . . . , 4 [26]. e average speed limits are considered for every route. All three highways that connect CT and CA are assumed to have 70 miles per hour speed limit, while the city roads are assumed to be 35 miles per hour, and the highway that connects Hamden and New Haven is assumed to be 55 miles per hour.
In addition, fuel costs (c u 1 ) and operational costs (c u 2 ) that include maintenance, repair, and depreciation are used to calculate the traveling costs for each route. From the American Automobile Association (AAA), we assume the fuel price to be $2.376 per gallon (AAA, 2019) and fuel consumption of 29 city/41 highway miles per gallon. We set the fuel consumption parameter, c u 1 , equal to $0.0813 per mile (i.e., $2.376/29). If a vehicle travels mostly on highways, then c u 1 is $0.0579 per mile (i.e., $2.376/41). To calculate the maintenance cost per mile, c u 2 , we use two data points: (i) the national, maintenance cost average for a new Sedan which is $242 per year (AAA, 2019); (ii) annual miles traveled on an average which is 13,476 [27]. Using (i) and (ii), we set c u 2 equal to $0.0179 (i.e., 242/13,476).
ere are three scenarios considered for each evaluation criteria of the objective functions z 1 , z 2 , z 3 , and z 5 : best, worst, and average. e best and worst scenarios correspond to the objective function values when the minimum and maximum duration are used for waiting times at the traffic signals, respectively. e average scenario corresponds to the mean of the objective function values when the waiting times are allowed to vary randomly between the minimum and maximum. For instance, for each objective function, we for u � 1 to b do input the sets K u , Δ u , δ u , V, y u 1 , and y u 2 and the output from Algorithm 4 (this input determines the sets Ψ u and P) Input the costs c 1 (u) and c 2 (u) Input the weights σ (1) h and σ (2) h for all h for y u 1 � 1 to p u 1 do (to calculate yield sign duration for a route) for y u 2 � 1 to p u 2 do μ u y 1 �� Rand(min μ u y 1 , max μ u y 1 ) μ u y 2 �� Rand(min μ u y 2 , max μ u y 2 ) end for end for for j � 1 to 4 do (to calculate stop sign duration for a route) w u j �� Rand(min w u j , max w u j ) end for for i � 1 to e do (to calculate traffic light duration for a route) t u i �� Rand(0, max β(i, u) ) end for for h � 1 to 5 do (to calculate z 6 and z 7 after receiving user preferences) input σ (1) h and σ (2) h calculate z h Sort z 6 and z 7 from min to max by using Quick sort algorithm display min z 1 −z 7 values with the corresponding specific route u choices end for end for ALGORITHM 5: Optimization of objective functions z 1 − z 7 . perform 1000 replications of simulation of the waiting times and collect the average statistics. e best, worst, and average scenarios are different for z 1 ; Algorithms 1-4 are used to calculate the total distance values for z 1 . is calculation has two phases: for instance, in the first phase, for every geodesic and linked chain circle combination, the minimum length of the geodesic (i.e., the minimum distance to travel from intersection to intersection) is determined; in the second phase, using the minimum length geodesics, the shortest path from O to D is identified. e best case value of z 1 is calculated using the shortest path geodesic (i.e., inner lane) distance, while the worst-case value is calculated using the longest geodesic distance (i.e., outer lane) on each route.
Since z 3 includes the cost incurred due to waiting at traffic signals, the best-and worst-case values of z 3 correspond to the shortest and longest duration of waiting times at the traffic signals. Similarly, the waiting times at the traffic signals are included in the total travel time; as a result, bestand worst-case values of z 5 correspond to the shortest and longest duration of waiting times. Note that there are no best-and worst-case values for z 4 as the number of traffic signals does not vary on a given route. Next, we analyze the computational results from the short-distance case.

Short-Distance Routing: A Case Study.
For demonstration, we choose two locations, a coffee shop in Hamden and the train station in New Haven, where Hamden and New Haven are two cities in the state of Connecticut. ere are three routes between the two places: (A) via I-91 South-using highway; (B) via Whitney Avenue-using local roads; (C) via Ridge Road-using local roads. On each route, there is a certain number of traffic lights, yield signs, and stop signs. Table 6 summarizes the computational results for the routes A, B, and C concerning criteria z 1 to z 6 . e best values are given in bold. We compare the results with Google Maps when the relevant information is available. We make the following observations: (1) When z 1 is considered, Route B is a good choice in all scenarios: best, worst, and average. e difference between the best and the worst cases is the highest for Route A which is 0.54 miles (11.38-10.84), about 5% per trip. In other words, if Route A is selected for traveling, using the innermost lanes can generate savings on the distance of about 5.
(2) When z 2 is considered, Route C is a good choice in all scenarios. ough the simulated average values for all routes are close to each other, Route C has the minimum. Moreover, Route C has the minimum value in the best-and worst-case scenarios. e waiting time information is not available from Google Maps. Note that the simulated average values obtained for z 2 has a significant impact on total travel time.
(3) Route B is the most cost effective when z 3 is considered. e simulated average value for Route B is $0.727. When compared to Routes A and C, Route B generates an average cost savings of 0.185 and 0.166, respectively. With respect to z 4 , Route A is the best as it has a minimum number of traffic signs. (4) When z 5 is considered, Route C is a good choice in all scenarios. If the innermost lanes (i.e., best-case scenario) are used to travel, both Routes A and C can save travel time when compared to Google Maps of the order of 1 minute (i.e., 6%). In the average scenario, Route C is a good choice; however, Google Maps performs better with a travel time of 17 minutes. e difference between our results and Google Maps is due to the variability in the waiting times that impact the total travel time. Further, we note that when the total distance is considered, Route A (highway) is the worst choice; however, it does well when travel time is considered due to the minimum number of traffic signals and increased speed limits. Overall, the above computational experiments for shortdistance routes highlight the complexity of decision-making in evaluating routes when multiple criteria are considered. With criteria z 1 to z 6 , a commuter can pick and choose from several travel options based on individual preference. e geodesic distances calculated based on the RMS approach are on par with the state-of-the-art application, Google Maps. By including traffic sign duration and geodesic distances, our approach enhances the Google Maps' results. Next, we analyze the computational results from long-distance problem instances.  Table 7 summarizes the computational results for the routes A, B, and C with respect to criteria z 1 to z 6 . e best values are given in bold. We compare the results with Google Maps when the relevant information is available. e following are the observations extracted from Table 7: (1) Route A is a good choice in the best-case scenario (using innermost lanes) when z 1 is considered. On the other hand, in the worst (using outermost lanes) and average (simulation) scenarios, Route B is a good choice. When compared to Google's best, Route A in the best-case scenario which generates savings of 24 miles (i.e., 3003 − 2978.88) per trip. ough the savings per trip is minimal, when a large fleet of trucks and hundreds of deliveries are considered, the total savings can be a considerable amount. Further, the difference between the best-and worst-case scenarios is at the maximum with Route A (3029.88 − 2978.88 � 51 miles).
is difference impacts the total travel time, cost of traveling an additional 51 miles, and all the costs associated with the transportation (e.g., driver, truck maintenance, and fuel consumption costs).
(2) Route C is the best choice in all scenarios when the total waiting time at traffic signs, z 2 , is considered.  the solution vectors and assist the decision maker(s) to examine the best options for all routes.
In summary, the geodesic distances calculated based on the RMS approach with traffic signals considered for the routing solution outperform the state-of-the-art application, Google Maps. ough our model does not consider the real-time traffic conditions and congestion effects, the results show that the geodesic distances calculated based on the RMS approach can be helpful in enhancing the accuracy of the results from Google Maps, see Figure 4. For instance, using the innermost lanes (best case) on Route A, one can generate distance savings of 25 miles per trip when compared to the best result from Google Maps. Google Maps displays a single distance for traveling from the origin to destination, while the nature of the traffic roads suggests a distance to be traveled between two values: minimum and maximum geodesic distances.
ough the savings may seem minimal, when the scope is increased from one trip to hundreds of trips per year and  Journal of Advanced Transportation from one truck to the entire fleet of trucks, the magnitude of the savings can be significant.

Conclusion and Future Research
Transportation is critical for both end consumers and the business world. One of the key components of the transportation is to identify the best route to travel that connects an origin and a destination. For instance, a route can be evaluated using multiple criteria such as travel cost and travel time. So, for some, the "best" route is cost effective, and for others, it is time sensitive. e approach to identify the routes, the assumptions, and the underlying calculations is important as they determine the accuracy of the routing solutions. In this research, we consider several real-world driving factors such as the time spent at traffic signs (e.g., yield signs and stop signs), speed limits, and the topology of the surface to develop realistic and accurate routing solutions. ough these factors increase the complexity of modeling, they provide the flexibility to evaluate the routing solutions from different perspectives: cost, distance, and time, to name a few.
In this research, we develop a nonlinear, multiobjective, mixed-integer mathematical model (MINLP) with five objectives that minimizes travel distance, travel time, travel cost, time spent on traffic signs, and the number of traffic signs to design and evaluate the routes between an origin and a destination where the waiting times associated with traffic lights, stop signs, and yield signs are stochastic. In addition, we present a weighted sum formulation (WSF) and a best vector set (BV) approach to assist in evaluating routing solutions. A set of algorithms are developed to solve the MINLP that are based on Riemannian manifold surfaces (RMS) to factor in the Earth's curvature when calculating distances. A modified version of the linked chain method (LCM) is used to determine geodesic distances that are important in calculating the distances for possible routes between an origin and a destination. e geodesic distances determined in the first four algorithms based on RMS are used in the MINLP model to ensure that the routing solutions are realistic and accurate. For instance, using the LCM, we generate inner-and outermost lanes for highways that provide a range of traveling distances (instead of one distance) which has not been addressed in the past research.
We demonstrate the complexity of our algorithms and apply them in real-life instances. Two different real-life routing cases, short and long distances, for computational experiments and discrete event simulations are used for solving the MINLP and WSF and determining BV. e first case includes three different short-distance routes, and the second includes three long-distance routes. In both cases, we show that the geodesic distances calculated based on the RMS approach are on par with the state-of-the-art application, Google Maps. For instance, using the innermost lanes (best case) can generate savings on the distance traveled when compared to the best result from Google Maps. When a large fleet of trucks and hundreds of deliveries are considered, the total savings can be significant. In addition, we show that the MINLP, WSF, and BV assist decision-making in evaluating routes when multiple criteria are considered. We believe that there are significant research opportunities for our work to be applied in the areas of traffic planning/management for routing decision strategies. e results obtained in this work highlight the impact of geodesic distance calculations and the role of traffic signs in evaluating the routes which cannot be overlooked. We show that the waiting times at the traffic signs can be significant and play a critical role in evaluating the routes for short-distance transportation problems. For instance, Route B which uses the local roads to connect a coffee shop and a train station minimizes the total travel distance. But, the same route is a bad choice when the travel time is considered since the route has many traffic signs that create travel interruptions. On the other hand, for long-distance transportation problems, geodesic distance calculations play a major role. For instance, Route A, an actual truckload route of a multinational stainless steel and specialty metal manufacturer located in the state of Connecticut that follows the interstate highway I-80, has the minimum travel distance of 2,978.88 miles based on the shortest geodesic. When compared with the longest geodesic distance of 3029.33 miles, Route A generates a savings of 50 miles per trip per truck (1.7%) which can be significant when the total number of trips and the fleet size are considered.
Our research has the following practical implications: (i) it is difficult to obtain the realistic and accurate routing solutions from the available mapping software when they do not consider the nature of the Earth's surface. Our approach that is based on the Riemannian manifold surface addresses this issue by considering the shortest path geodesics. (ii) Given the magnitude of the transportation spending, i.e., 0.91 trillion dollars in 2014, and the pervasiveness of the distance calculations, our approach can be applied in many different areas such as the supply chain (e.g., network design), transportation (e.g., vehicle routing), renewable energy (e.g., market distribution), and alternative technology (e.g., location-allocation of charging stations for electric vehicles), to name a few. (iii) Furthermore, our multicriteria mathematical model where we include the waiting times at the traffic signs to evaluate the routes based on distance, time, and cost addresses the lack of consideration of the traffic signs in the state-of-the-art applications (e.g., Google Maps).

Data Availability
A set of data is collected by using the distance calculator of Google Maps and locations readily available in Google Maps; therefore, the numerical results attained in this work as a part of the tables can be attained by other researchers by using Google Maps. e data set and the corresponding analysis can be shared by Dr. Emre Tokgoz via e-mail (Emre.Tokgoz@qu.edu) if requested.

Conflicts of Interest
e authors declare that they have no conflicts of interest.