# Advancing Transportation Routing Decisions Using Riemannian Manifold Surfaces

Authors: Emre Tokgöz; Iddrisu Awudu; Saravanan Kuppusamy
Journal: Journal of Advanced Transportation
 (2020)

---

## Abstract

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. Though 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. 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 traffic signs to design and evaluate the routes where the waiting times associated with traffic 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 significant cost savings.


---

## Body

## 1. 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. The distance between locations is often used as a surrogate measure for the transportation cost. It is therefore vital to understand how distances are calculated.The 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., ageodesic. The 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. The 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.


### 1.1. 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. There 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. Figure1 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. The dashed lines in Figure 1 represent the inner- and outermost lanes of the road. The 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 RMS-based approach (the details are provided later in the paper). The 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. This 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.Figure 1
An illustration of RMS.Table 1
Summary for routes A, B, and C.
RoutesDistance-based outermost lanes (miles)Distance-based innermost lanes (miles)Savings (miles)A11.3810.840.53B7.227.210.01C7.3667.3670.001Next, 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. Table2 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.Table 2
Evaluation summary for routes A, B, and C.
RoutesDistance (RMS) (miles)Waiting time at traffic signs (mins)Travel costTotal number of traffic signsTraveling time (mins)A10.844.1$0.8681815.93B7.215.65$0.7152918.01C7.373.3$0.7292115.93For 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. The 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.The rest of this paper is organized as follows: Section2 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.


## 1.1. 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. There 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. Figure1 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. The dashed lines in Figure 1 represent the inner- and outermost lanes of the road. The 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 RMS-based approach (the details are provided later in the paper). The 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. This 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.Figure 1
An illustration of RMS.Table 1
Summary for routes A, B, and C.
RoutesDistance-based outermost lanes (miles)Distance-based innermost lanes (miles)Savings (miles)A11.3810.840.53B7.227.210.01C7.3667.3670.001Next, 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. Table2 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.Table 2
Evaluation summary for routes A, B, and C.
RoutesDistance (RMS) (miles)Waiting time at traffic signs (mins)Travel costTotal number of traffic signsTraveling time (mins)A10.844.1$0.8681815.93B7.215.65$0.7152918.01C7.373.3$0.7292115.93For 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. The 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.The rest of this paper is organized as follows: Section2 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.


## 2. 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.


### 2.1. 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. Therefore, 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 lp-norm, which has been proven to be an accurate distance calculation function, to the l2b-norm which uses different parameter estimates for the distance prediction. This 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.


### 2.2. 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. Though 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. Therefore, 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 algorithm-based 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.


### 2.3. Research Contributions of This 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. The authors conclude that the skewed transportation networks in large-scale 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. This 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, mixed-integer 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).


## 2.1. 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. Therefore, 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 lp-norm, which has been proven to be an accurate distance calculation function, to the l2b-norm which uses different parameter estimates for the distance prediction. This 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.


## 2.2. 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. Though 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. Therefore, 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 algorithm-based 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.


## 2.3. Research Contributions of This 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. The authors conclude that the skewed transportation networks in large-scale 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. This 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, mixed-integer 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).


## 3. Riemannian Manifold Surfaces and Linked Chain Method

Carl Friedrich Gauss provedTheorema Egregium (meaning remarkable theorem in Latin) in 1828 by establishing an important property of surfaces. The 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. The 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 positive-definite quadratic forms on the tangent bundle when the length of a continuously differentiable curve γ:a,b⟶M is defined by the following equation:(1)lγ=∫a,bdγtdtdt.The 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γ when γ 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.The 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. This change in the curvature is important when calculating shortest path distances on Earth’s surface. This surface is a natural manifold, and the pathways on this manifold are natural geodesics.The 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. The 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 inner- and 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. The 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.


## 4. 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. The set of notations used in the development of algorithms is given in Table 3.Table 3
Notations.
MManifold surfaceOOrigin onMDDestination onMM0Predetermined neighborhood ofO−D that contains both O and DmaxM0Maximum distance from the origin,O, to the boundary of M0rRadius of LCM circlesr0An initially chosen constant radius determined by the decision maker (s)γirTheith geodesic identified with radius rCγirThe LCM circle that contains the geodesicγirnrNumber of LCM circles with radiusr that covers M0ϵA constant that is used to increment the radius of a LCM circleGSet of geodesicsPSet of all vectors with each vector containing a geodesic, length of the geodesic, and the LCM circle that contains the geodesicNSet of linked list vectors formed by the geodesicsWe assume geodesics are determined based on the pathway intersections and illustrate this assumption by using a T-intersection. For instance, a T-intersection, see Figure2, 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. The 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.,(2)minz1=minlγ=∫O,Ddγtdtdt.Figure 2
An illustration of geodesics using the T-intersection.


### 4.1. Subroutine 1

Algorithm1 shows the pseudocode for the first subroutine. First, a surrounding area M0 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., γir’s, they are added to a set G. This 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.Algorithm 1: Geodesic search pseudocode.
Input: (a) coordinates of O, D, and M0; (b) GIS map of the roads within M0 on M.Run: (a) G=∅; (b) r=r0, where r0 = initial radius; (c) ϵ=0.001.(1)
forr=r0 to maxM0do(2)
CoverM0 with LCM-linked circles with radius r by tracking all possible roads between O and D(3)
fori=1 to nrdo(4)
Assignγir to be the ith geodesic covered for radius r(5)
AssignCγir to be the circle that covers γir(6)
ifγir∉Gthen(7)
G⟵G∪γir(8)
else(9)
r=r+ϵ(10)
end if(11)
end for(12)
end forAlgorithm 2: Geodesic identifier and distance calculation pseudocode.
Run: Algorithm 1Input: (a) set G from Algorithm 1; (b) Gnew=∅=P.(1)
forr=r0 to maxM0do(2)
fori=1 to nrdo(3)
forγir in Gdo(4)
Calculatelγir(5)
Definep^γir=γir,lγir,Cγir(6)
Gnew⟵Gnew∪γir(7)
G⟵G−γir(8)
P⟵P∪p^γir(9)
end for(10)
end for(11)
end for(12)
AssignG=GnewFigure3 illustrates the steps involved in Algorithm 1. For illustration purposes, a focus area M0 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 M0 is covered with the circles of initial radius r0. For the purposes of clarity, the fully covered M0 in 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γ1r0, Cγ2r0, Cγ3r0, and Cγ4r0 to geodesics γ1r0, γ2r0, γ3r0, and γ4r0, respectively. Finally, the above geodesics are added to set G. The above steps are repeated with a new radius, r0=r0+ϵ. If geodesics are identified with the circles of new radii, they will be added to G, otherwise, not.Figure 3
An illustration of Algorithm1. (a) Focus area with r = r0. (b) Focus area with r = r1. (c) Focus area with multiple radii.
(a)(b)(c)Figure3(b) illustrates the geodesics that are identified when the radius is equal to r1. In this figure, the circles that identify geodesics are represented by the dashed lines. We assign circles Cγ1r1, Cγ2r1, Cγ3r1, and Cγ4r1 to geodesics γ1r1, γ2r1, γ3r1, and γ4r1, respectively. We add the above geodesics to set G. Figure 3(c) is an illustration of the focus area when all geodesics are identified. In this figure, 22 geodesics are identified, which are numbered from 1 to 22 by 22 different circles which are in G.


### 4.2. Subroutine 2

Algorithm2 shows the pseudocode for the second subroutine. We use this subroutine to calculate the length of geodesics, γ, that are in G. Initially, two empty sets are created, Gnew and P. First, in set G, we select a radius r and examine the geodesics, γir, identified with that radius. Second, the length, lγir, of every geodesic is calculated using the following function:(3)lγir=∫Cγirdγirtdtdt.Third, we definep^γir to be a vector with three elements: the geodesic γir, the length of the geodesic lγir, and the circle Cγir that covers the geodesic. Subsequently, we add the geodesic γir to set Gnew and delete the same geodesic from set G. Finally, p^γir is added to set P. The 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.


### 4.3. Subroutine 3

In topology, geodesics are connected and form a path if circles that contain them are overlapping. This subroutine searches for overlapping circles to determine the paths within the focus area,M0. 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. Then, two elements p^γir and p^γjr of set P are chosen such that i≠j. If circles Cir and Cjr that contain γir and γjr are overlapping, the link is added to set N; otherwise, a different pair of p^ is chosen. This 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.Algorithm 3: Geodesic path construction pseudocode.
Run: Algorithm 2InputN=∅(1)
fori=1 to nrdo(2)
whilep^γir∈Pdo(3)
forj=1 to nrdo(4)
ifi≠j and Cγir∩Cγjr≠∅then(5)
Linkγir and γjr and add the link to N(6)
Continue updatingN with links for p^γir∈P(7)
end if(8)
end for(9)
end while(10)
end for


### 4.4. Subroutine 4

Algorithm4 shows the pseudocode for the fourth subroutine. First, P is examined to identify the circles that contain origin, O, and destination, D. Once the circles are identified, they are split into nonoverlapping circles and geodesics. This is done to improve the accuracy of calculation of the distance between O and D. For instance, if Co is the circle that contains O and γo is the geodesic in Co, then Co is split into Coright and Coleft and γo is split into γoright and γoleft. Subsequently, the corresponding vectors, p^γoright and p^γoleft, are added to P. Similarly, Cd, the circle that contains D, is split into Cdright and Cdleft and γd is split into γdright and γdleft and the corresponding vectors, p^γdright and p^γdleft, are added to P. We assume h number of pathways exist between O and D, and the set Nh, a subset of N, contains all those pathways. The length, lγh, of all pathways is calculated to determine all possible distances from O to D. Finally, the path with the least distance, lmin, is selected to be the shortest path between O and D.Algorithm 4: Path finder and distance minimizer pseudocode.
Run Algorithm 3(1)
forp^γir∈Pdo(2)
ifO∈Cγirthen(3)
Assignγo:=γir and Co:=Cγir(4)
SplitCo into Coright and Coleft and γo into γoright and γoleft(5)
Calculatelγoleft and lγoright(6)
Definep^γoright=γoright,lγoright,Cγoright and p^γoleft=γoleft,lγoleft,Cγoleft(7)
P⟵P∪p^γoright,p^γoleft(8)
end if(9)
ifD∈Cγirthen(10)
Repeat steps 3–6 by replacingO with D(11)
P⟵P∪p^γdright,p^γdleft(12)
end if(13)
ifp^γoright,p^γoleft∪p^γdright,p^γdleft∈Nh⊂Nthen(14)
Find the minimum and maximum oflγir calculated in Algorithm 2(15)
Calculate the total length of the geodesics for allNh using step 14(16)
lmin=minNh and lmax=maxNh are the minimum and maximum length geodesics for all Nh that connect O and D(17)
end if(18)
end forNote that steps 14–16 in Algorithm4 are used to calculate minimum and maximum geodesic distances lγir for all γ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.
There 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. Therefore, f is an onto map; there exists a 1-1 and onto map that projects C to G. The inverse map f−1:G⟶C maps each γir in G to Cγir in C.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. Algorithm2 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 path-connected predetermined region. The set Nh 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, Algorithms1–4 are not applicable.


## 4.1. Subroutine 1

Algorithm1 shows the pseudocode for the first subroutine. First, a surrounding area M0 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., γir’s, they are added to a set G. This 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.Algorithm 1: Geodesic search pseudocode.
Input: (a) coordinates of O, D, and M0; (b) GIS map of the roads within M0 on M.Run: (a) G=∅; (b) r=r0, where r0 = initial radius; (c) ϵ=0.001.(1)
forr=r0 to maxM0do(2)
CoverM0 with LCM-linked circles with radius r by tracking all possible roads between O and D(3)
fori=1 to nrdo(4)
Assignγir to be the ith geodesic covered for radius r(5)
AssignCγir to be the circle that covers γir(6)
ifγir∉Gthen(7)
G⟵G∪γir(8)
else(9)
r=r+ϵ(10)
end if(11)
end for(12)
end forAlgorithm 2: Geodesic identifier and distance calculation pseudocode.
Run: Algorithm 1Input: (a) set G from Algorithm 1; (b) Gnew=∅=P.(1)
forr=r0 to maxM0do(2)
fori=1 to nrdo(3)
forγir in Gdo(4)
Calculatelγir(5)
Definep^γir=γir,lγir,Cγir(6)
Gnew⟵Gnew∪γir(7)
G⟵G−γir(8)
P⟵P∪p^γir(9)
end for(10)
end for(11)
end for(12)
AssignG=GnewFigure3 illustrates the steps involved in Algorithm 1. For illustration purposes, a focus area M0 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 M0 is covered with the circles of initial radius r0. For the purposes of clarity, the fully covered M0 in 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γ1r0, Cγ2r0, Cγ3r0, and Cγ4r0 to geodesics γ1r0, γ2r0, γ3r0, and γ4r0, respectively. Finally, the above geodesics are added to set G. The above steps are repeated with a new radius, r0=r0+ϵ. If geodesics are identified with the circles of new radii, they will be added to G, otherwise, not.Figure 3
An illustration of Algorithm1. (a) Focus area with r = r0. (b) Focus area with r = r1. (c) Focus area with multiple radii.
(a)(b)(c)Figure3(b) illustrates the geodesics that are identified when the radius is equal to r1. In this figure, the circles that identify geodesics are represented by the dashed lines. We assign circles Cγ1r1, Cγ2r1, Cγ3r1, and Cγ4r1 to geodesics γ1r1, γ2r1, γ3r1, and γ4r1, respectively. We add the above geodesics to set G. Figure 3(c) is an illustration of the focus area when all geodesics are identified. In this figure, 22 geodesics are identified, which are numbered from 1 to 22 by 22 different circles which are in G.


## 4.2. Subroutine 2

Algorithm2 shows the pseudocode for the second subroutine. We use this subroutine to calculate the length of geodesics, γ, that are in G. Initially, two empty sets are created, Gnew and P. First, in set G, we select a radius r and examine the geodesics, γir, identified with that radius. Second, the length, lγir, of every geodesic is calculated using the following function:(3)lγir=∫Cγirdγirtdtdt.Third, we definep^γir to be a vector with three elements: the geodesic γir, the length of the geodesic lγir, and the circle Cγir that covers the geodesic. Subsequently, we add the geodesic γir to set Gnew and delete the same geodesic from set G. Finally, p^γir is added to set P. The 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.


## 4.3. Subroutine 3

In topology, geodesics are connected and form a path if circles that contain them are overlapping. This subroutine searches for overlapping circles to determine the paths within the focus area,M0. 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. Then, two elements p^γir and p^γjr of set P are chosen such that i≠j. If circles Cir and Cjr that contain γir and γjr are overlapping, the link is added to set N; otherwise, a different pair of p^ is chosen. This 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.Algorithm 3: Geodesic path construction pseudocode.
Run: Algorithm 2InputN=∅(1)
fori=1 to nrdo(2)
whilep^γir∈Pdo(3)
forj=1 to nrdo(4)
ifi≠j and Cγir∩Cγjr≠∅then(5)
Linkγir and γjr and add the link to N(6)
Continue updatingN with links for p^γir∈P(7)
end if(8)
end for(9)
end while(10)
end for


## 4.4. Subroutine 4

Algorithm4 shows the pseudocode for the fourth subroutine. First, P is examined to identify the circles that contain origin, O, and destination, D. Once the circles are identified, they are split into nonoverlapping circles and geodesics. This is done to improve the accuracy of calculation of the distance between O and D. For instance, if Co is the circle that contains O and γo is the geodesic in Co, then Co is split into Coright and Coleft and γo is split into γoright and γoleft. Subsequently, the corresponding vectors, p^γoright and p^γoleft, are added to P. Similarly, Cd, the circle that contains D, is split into Cdright and Cdleft and γd is split into γdright and γdleft and the corresponding vectors, p^γdright and p^γdleft, are added to P. We assume h number of pathways exist between O and D, and the set Nh, a subset of N, contains all those pathways. The length, lγh, of all pathways is calculated to determine all possible distances from O to D. Finally, the path with the least distance, lmin, is selected to be the shortest path between O and D.Algorithm 4: Path finder and distance minimizer pseudocode.
Run Algorithm 3(1)
forp^γir∈Pdo(2)
ifO∈Cγirthen(3)
Assignγo:=γir and Co:=Cγir(4)
SplitCo into Coright and Coleft and γo into γoright and γoleft(5)
Calculatelγoleft and lγoright(6)
Definep^γoright=γoright,lγoright,Cγoright and p^γoleft=γoleft,lγoleft,Cγoleft(7)
P⟵P∪p^γoright,p^γoleft(8)
end if(9)
ifD∈Cγirthen(10)
Repeat steps 3–6 by replacingO with D(11)
P⟵P∪p^γdright,p^γdleft(12)
end if(13)
ifp^γoright,p^γoleft∪p^γdright,p^γdleft∈Nh⊂Nthen(14)
Find the minimum and maximum oflγir calculated in Algorithm 2(15)
Calculate the total length of the geodesics for allNh using step 14(16)
lmin=minNh and lmax=maxNh are the minimum and maximum length geodesics for all Nh that connect O and D(17)
end if(18)
end forNote that steps 14–16 in Algorithm4 are used to calculate minimum and maximum geodesic distances lγir for all γ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.
There 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. Therefore, f is an onto map; there exists a 1-1 and onto map that projects C to G. The inverse map f−1:G⟶C maps each γir in G to Cγir in C.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. Algorithm2 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 path-connected predetermined region. The set Nh 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, Algorithms1–4 are not applicable.


## 5. 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 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. The set of notations used in the mathematical model is given in Table 4.Table 4
A summary of model notations.
bNumber of routes that connects the origin (O) and the destination (D)KuSet of 1-way, 2-way, 3-way, and 4-way stop signs on routeuKjuSet ofj‐way stop signs on route uηjuNumber ofj‐way stop signs on route uwjWait time at aj‐way stop signI0Number of traffic light categoriesΔuSet of traffic lights on routeuΔiuSet of theith category traffic lights on route uδiuNumber of theith category of traffic lights on route uδuTotal number of traffic lights of all categories on routeutiuWait time at anith category traffic light on route uy1uNumber of yield signs to enter highway on routeuy2uNumber of yield signs to enter city traffic on routeuμykuWait time at theykth sign on route uVSet ofq-types of speed limits from O to DγsuGeodesics on route u that connects O and DΨuSet of all possible geodesics on routeuc1uFixed per-mile fuel cost on routeuc2uVariable per-mile operational cost on routeuzhObjective function with criteriah, where h = 1 (distance), =2 (wait time), =3 (cost), =4 (number of traffic signs), =5 (travel time), and =6 (weighted sum)σh1Assigned weight for criteriah


### 5.1. Assumptions and Parameters

(i)
Stop Signs. There are four types of stop signs used for mathematical formulation: 1-way, 2-way, 3-way, and 4-way. Let ηju be the number of stop signs for the j‐way stop on route u for all 1≤u≤b, where b is the total number of routes that connect O to D. The 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, wj, a stochastic parameter, varies between waiting times minwj and maxwj. The set of all stop signs is(4)Ku=K1u,K2u,K3u,K4u,where(5)Kju=kj−1u+1,kj−1u+2,…,kju,represents the set ofj−way stop signs kj−1u 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 minwj=3 and maxwj=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(6)Δiu=δi−1u+1,δi−1u+2,…,δiu,represent the category-i traffic lights based on βiu 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, Δ1100 represents the 1st category of traffic lights that has a maximum duration of 25 seconds on the 100th route. The set of traffic lights for b different paths from O to D takes the form(7)Δu=Δ1u,Δ2u,…,Δeu,for all1≤u≤b and e∈I0 where I0 is the number of traffic light categories. The number of traffic lights on the uth pathway in category i is represented by Δiu, and the total number of traffic lights on each path u from O to D is(8)δu=∑i=1eΔiu,for all1≤u≤b. For the example given above, if the set Δ1100 has 5 elements, then Δ1u=5 indicates that there are 5 traffic lights with 25 seconds of maximum wait time on route 100. Let tiu represent the stochastic parameter that varies between 0 and βiu for waiting at the category-i traffic light on route u.(iii)
Yield Signs. There 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 y1u represent the number of yield signs to enter the highway and y2u represent the number of yield signs to enter the city traffic on route u for all 1≤u≤b. Let yu=y1u+y2u for each u, and let μyku be the stochastic parameter representing the wait time at yku 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. The distance calculations depend on speed limit, geodesic distance, and traffic sign durations. The 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(9)V=v1,v2,…,vq,where there areq 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 γs from Algorithm 2 (subscript ir is modified as s) and the corresponding length lγs, the distance traveled on the LCM circle Cγs can be calculated as follows: lγs/vl, where vl∈V is the speed limit value on the road covered by Cγs. Later, we use this formula in the MIP formulations to determine the travel time from O to D.


### 5.2. Objective Functions and Constraints

In this section, we introduce five objective functions and develop three different formulations. Letγsu represent the geodesic s on route u that connects O and D and Ψu=γsus∈I1 is a set of all possible geodesics between O and D with the index set I1. Let c1u and c2u be the fixed per-mile fuel cost and per-mile operational cost, respectively. The constrained, mixed-integer, nonlinear program formulation for the b-path problem to travel from O to D that includes z1 can be formulated as follows:(10)total distanceminz1=∫O,Ddγsdsd,(11)total waiting timeminz2=∑i=1etiuΔiu+∑j=14wjuηju+∑k=12ykuμyku,(12)total costminz3=c1u∫O,Ddγsu+c2uz2u,(13)total number of traffic signsminz4=∑i=1eΔiu+∑j=14ηju+∑k=12yku,(14)total travel timeminz5=∑slγsuvlu,subject to(15)minμy1u≤μy1u≤maxμy1u,∀y1=1,…,y1uand1≤u≤b,(16)minμy2u≤μy2u≤maxμy2u,∀y2=1,…,y2uand1≤u≤b,(17)minwju≤wju≤maxwju,∀j=1,2,3,4;kj−1u∈Kuand1≤u≤b,(18)0≤tiu≤maxβiu,∀i∈I0=1,2,…,eand1≤u≤b,(19)γsu∈Ψu,∀s∈I1and1≤u≤b,(20)Zh=zhuu=1b,∀h=1,2,…,5,(21)vlu∈V,∀l∈I2and1≤u≤b.Recall that algorithms from Section4 are used to determine the shortest geodesic path, and the corresponding objective function z1 minimizes the travel distance. We reintroduce z1 in (10) for clarity. In this formulation, u is an integer variable that denotes the route. In addition, γt in (10) is a nonlinear, continuous variable that is used to calculate the geodesic distance on nonlinear curves. The 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 (11) minimizes the total waiting time on traffic signals while traveling from O to D. Equation (12) minimizes the cost incurred (in dollars) to travel from O to D. Equation (13) minimizes the total number of traffic signs, and equation (14) minimizes the total travel time. Equations (15) and (16) determine the waiting times at the city and highway yield signs, respectively. Equation (17) determines the waiting time at stop signs, and equation (18) provides the waiting times at the traffic lights. The set element specifications are given in (19)–(21). Equations (10)–(14) determine the minimum values of zhu for all h and u.


#### 5.2.1. Weighted Sum Formulation (WSF)

The weighted sum formulation (WSF) allows the decision maker to choose a route based on the weights assigned to the objective functionsz1 through z5. Function f. introduced below normalizes the objective function values z1–z5 for each route, and it is used as a part of minimizing the weighted sum. This option of formulation requires weighted sum formulation of equations (15)–(19) by introducing the weight function σh:A1⟶ℝ+ with the corresponding objective function:(22)minz6=∑h=15σh1fzhu,where A1 is the finite index set and fzhu are the normalized values of the objective functions zhu calculated for each route u by using the formula:(23)fzhu=zhu∑j=1bzhj,for each h=1,…,5 and each u such that 1≤u≤b. The values of fzhu are determined for all h and u that forms a set; this set of values is used for determining the minimum value of z6. The weights σh1 are known constants that are assigned by the decision maker, and the sum of the weights must be equal to 100%.


#### 5.2.2. 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 eachzi,∀i=1,…,5. Since the optimal route can vary with zi, the corresponding vector set solution can be different. In this case,(24)minz7u=z1u,…,minziu,minzju,…,z5u,where i,j=1,…,5, u=1,…,b, b is the total number of the routes, and the best vector set BV=minz7u which is a subset of all solutions. For instance, assume b = 4. Let route 1 (i.e., u = 1) be optimal for z1 and z2, route 2 (u = 2) be optimal for z3 and z5, route 3 (u = 3) be optimal for z3, z4, and z5, and route 4 (u = 4) be not optimal for any zi. The best vector (BV) solution is given as follows:(25)BV=minz11,minz21,z31,z41,z51,z12,z22,minz32,z42,minz52,z13,z23,minz33,minz43,minz53.Note in BV that there is no minimum objective function value for route 4 as it is not optimal for anyzi. The 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.


### 5.3. An Algorithmic Solution

In this section, we develop an algorithmic solution for the mathematical models presented in Section5.2. Recall that equations (10), (12), and (14) are nonlinear since the distance calculations of geodesics depend on nonlinear curves. The numbers of traffic controls such as stop signs, yield signs, and traffic lights are known constants for each route. Therefore, the objective function in equation (13) is a constant for a route. Algorithm 5 uses the outputs of Algorithms 1–4 from Section 4.Algorithm 5: Optimization of objective functions z1−z7.
foru = 1 to bdoinput the setsKu, Δu, δu, V, y1u,andy2u and the output from Algorithm 4 (this input determines the sets Ψu and P)Input the costsc1u and c2uInput the weightsσh1 and σh2 for all hfory1u=1 to p1udo (to calculate yield sign duration for a route)fory2u=1 to p2udoμy1u==Randminμy1u,maxμy1uμy2u==Randminμy2u,maxμy2uend forend forforj = 1 to 4 do (to calculate stop sign duration for a route)wju==Randminwju,maxwjuend forfori = 1 to e do (to calculate traffic light duration for a route)tiu==Rand0,maxβi,uend forforh = 1 to 5 do (to calculatez6andz7after receiving user preferences)inputσh1 and σh2calculatezhSortz6 and z7 from min to max by using Quick sort algorithmdisplay minz1−z7 values with the corresponding specific route u choicesend forend forIn 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. These constraint values are then used for calculating z2−z5. Values of σh are the inputs to the algorithm for calculating the weighted sum value z6. The values attained from z1−z5 are then used for calculating the best vector set solution z7; the decision maker can choose a route from a variety of routes by using z7 that depends on the minimum values attained for z1−z5. A real-life application of Algorithm 5 with specific entries for minimum and maximum values is given in Section 6.For instance, the setP which is an output of Algorithm 2 is a critical input to Algorithm 5. The set P consists of p^γs vectors which has the following information: (i) the geodesic γs, (ii) the length of the geodesic lγs, and (iii) the LCM circle Cγs that has γs. The stochastic parameters are calculated by using a uniform distribution.Note that running Algorithm5 for sufficient number of times and averaging the outcomes will yield the desired objective function values. The highest and lowest values of the parameters presented in constraints (15)–(18) determine the minimum and maximum values for the objective function z2. The locations and numbers of traffic signs along with the speed limits and maximum durations of traffic lights are assumed to be known.


## 5.1. Assumptions and Parameters

(i)
Stop Signs. There are four types of stop signs used for mathematical formulation: 1-way, 2-way, 3-way, and 4-way. Let ηju be the number of stop signs for the j‐way stop on route u for all 1≤u≤b, where b is the total number of routes that connect O to D. The 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, wj, a stochastic parameter, varies between waiting times minwj and maxwj. The set of all stop signs is(4)Ku=K1u,K2u,K3u,K4u,where(5)Kju=kj−1u+1,kj−1u+2,…,kju,represents the set ofj−way stop signs kj−1u 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 minwj=3 and maxwj=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(6)Δiu=δi−1u+1,δi−1u+2,…,δiu,represent the category-i traffic lights based on βiu 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, Δ1100 represents the 1st category of traffic lights that has a maximum duration of 25 seconds on the 100th route. The set of traffic lights for b different paths from O to D takes the form(7)Δu=Δ1u,Δ2u,…,Δeu,for all1≤u≤b and e∈I0 where I0 is the number of traffic light categories. The number of traffic lights on the uth pathway in category i is represented by Δiu, and the total number of traffic lights on each path u from O to D is(8)δu=∑i=1eΔiu,for all1≤u≤b. For the example given above, if the set Δ1100 has 5 elements, then Δ1u=5 indicates that there are 5 traffic lights with 25 seconds of maximum wait time on route 100. Let tiu represent the stochastic parameter that varies between 0 and βiu for waiting at the category-i traffic light on route u.(iii)
Yield Signs. There 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 y1u represent the number of yield signs to enter the highway and y2u represent the number of yield signs to enter the city traffic on route u for all 1≤u≤b. Let yu=y1u+y2u for each u, and let μyku be the stochastic parameter representing the wait time at yku 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. The distance calculations depend on speed limit, geodesic distance, and traffic sign durations. The 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(9)V=v1,v2,…,vq,where there areq 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 γs from Algorithm 2 (subscript ir is modified as s) and the corresponding length lγs, the distance traveled on the LCM circle Cγs can be calculated as follows: lγs/vl, where vl∈V is the speed limit value on the road covered by Cγs. Later, we use this formula in the MIP formulations to determine the travel time from O to D.


## 5.2. Objective Functions and Constraints

In this section, we introduce five objective functions and develop three different formulations. Letγsu represent the geodesic s on route u that connects O and D and Ψu=γsus∈I1 is a set of all possible geodesics between O and D with the index set I1. Let c1u and c2u be the fixed per-mile fuel cost and per-mile operational cost, respectively. The constrained, mixed-integer, nonlinear program formulation for the b-path problem to travel from O to D that includes z1 can be formulated as follows:(10)total distanceminz1=∫O,Ddγsdsd,(11)total waiting timeminz2=∑i=1etiuΔiu+∑j=14wjuηju+∑k=12ykuμyku,(12)total costminz3=c1u∫O,Ddγsu+c2uz2u,(13)total number of traffic signsminz4=∑i=1eΔiu+∑j=14ηju+∑k=12yku,(14)total travel timeminz5=∑slγsuvlu,subject to(15)minμy1u≤μy1u≤maxμy1u,∀y1=1,…,y1uand1≤u≤b,(16)minμy2u≤μy2u≤maxμy2u,∀y2=1,…,y2uand1≤u≤b,(17)minwju≤wju≤maxwju,∀j=1,2,3,4;kj−1u∈Kuand1≤u≤b,(18)0≤tiu≤maxβiu,∀i∈I0=1,2,…,eand1≤u≤b,(19)γsu∈Ψu,∀s∈I1and1≤u≤b,(20)Zh=zhuu=1b,∀h=1,2,…,5,(21)vlu∈V,∀l∈I2and1≤u≤b.Recall that algorithms from Section4 are used to determine the shortest geodesic path, and the corresponding objective function z1 minimizes the travel distance. We reintroduce z1 in (10) for clarity. In this formulation, u is an integer variable that denotes the route. In addition, γt in (10) is a nonlinear, continuous variable that is used to calculate the geodesic distance on nonlinear curves. The 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 (11) minimizes the total waiting time on traffic signals while traveling from O to D. Equation (12) minimizes the cost incurred (in dollars) to travel from O to D. Equation (13) minimizes the total number of traffic signs, and equation (14) minimizes the total travel time. Equations (15) and (16) determine the waiting times at the city and highway yield signs, respectively. Equation (17) determines the waiting time at stop signs, and equation (18) provides the waiting times at the traffic lights. The set element specifications are given in (19)–(21). Equations (10)–(14) determine the minimum values of zhu for all h and u.


### 5.2.1. Weighted Sum Formulation (WSF)

The weighted sum formulation (WSF) allows the decision maker to choose a route based on the weights assigned to the objective functionsz1 through z5. Function f. introduced below normalizes the objective function values z1–z5 for each route, and it is used as a part of minimizing the weighted sum. This option of formulation requires weighted sum formulation of equations (15)–(19) by introducing the weight function σh:A1⟶ℝ+ with the corresponding objective function:(22)minz6=∑h=15σh1fzhu,where A1 is the finite index set and fzhu are the normalized values of the objective functions zhu calculated for each route u by using the formula:(23)fzhu=zhu∑j=1bzhj,for each h=1,…,5 and each u such that 1≤u≤b. The values of fzhu are determined for all h and u that forms a set; this set of values is used for determining the minimum value of z6. The weights σh1 are known constants that are assigned by the decision maker, and the sum of the weights must be equal to 100%.


### 5.2.2. 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 eachzi,∀i=1,…,5. Since the optimal route can vary with zi, the corresponding vector set solution can be different. In this case,(24)minz7u=z1u,…,minziu,minzju,…,z5u,where i,j=1,…,5, u=1,…,b, b is the total number of the routes, and the best vector set BV=minz7u which is a subset of all solutions. For instance, assume b = 4. Let route 1 (i.e., u = 1) be optimal for z1 and z2, route 2 (u = 2) be optimal for z3 and z5, route 3 (u = 3) be optimal for z3, z4, and z5, and route 4 (u = 4) be not optimal for any zi. The best vector (BV) solution is given as follows:(25)BV=minz11,minz21,z31,z41,z51,z12,z22,minz32,z42,minz52,z13,z23,minz33,minz43,minz53.Note in BV that there is no minimum objective function value for route 4 as it is not optimal for anyzi. The 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.


## 5.2.1. Weighted Sum Formulation (WSF)

The weighted sum formulation (WSF) allows the decision maker to choose a route based on the weights assigned to the objective functionsz1 through z5. Function f. introduced below normalizes the objective function values z1–z5 for each route, and it is used as a part of minimizing the weighted sum. This option of formulation requires weighted sum formulation of equations (15)–(19) by introducing the weight function σh:A1⟶ℝ+ with the corresponding objective function:(22)minz6=∑h=15σh1fzhu,where A1 is the finite index set and fzhu are the normalized values of the objective functions zhu calculated for each route u by using the formula:(23)fzhu=zhu∑j=1bzhj,for each h=1,…,5 and each u such that 1≤u≤b. The values of fzhu are determined for all h and u that forms a set; this set of values is used for determining the minimum value of z6. The weights σh1 are known constants that are assigned by the decision maker, and the sum of the weights must be equal to 100%.


## 5.2.2. 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 eachzi,∀i=1,…,5. Since the optimal route can vary with zi, the corresponding vector set solution can be different. In this case,(24)minz7u=z1u,…,minziu,minzju,…,z5u,where i,j=1,…,5, u=1,…,b, b is the total number of the routes, and the best vector set BV=minz7u which is a subset of all solutions. For instance, assume b = 4. Let route 1 (i.e., u = 1) be optimal for z1 and z2, route 2 (u = 2) be optimal for z3 and z5, route 3 (u = 3) be optimal for z3, z4, and z5, and route 4 (u = 4) be not optimal for any zi. The best vector (BV) solution is given as follows:(25)BV=minz11,minz21,z31,z41,z51,z12,z22,minz32,z42,minz52,z13,z23,minz33,minz43,minz53.Note in BV that there is no minimum objective function value for route 4 as it is not optimal for anyzi. The 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.


## 5.3. An Algorithmic Solution

In this section, we develop an algorithmic solution for the mathematical models presented in Section5.2. Recall that equations (10), (12), and (14) are nonlinear since the distance calculations of geodesics depend on nonlinear curves. The numbers of traffic controls such as stop signs, yield signs, and traffic lights are known constants for each route. Therefore, the objective function in equation (13) is a constant for a route. Algorithm 5 uses the outputs of Algorithms 1–4 from Section 4.Algorithm 5: Optimization of objective functions z1−z7.
foru = 1 to bdoinput the setsKu, Δu, δu, V, y1u,andy2u and the output from Algorithm 4 (this input determines the sets Ψu and P)Input the costsc1u and c2uInput the weightsσh1 and σh2 for all hfory1u=1 to p1udo (to calculate yield sign duration for a route)fory2u=1 to p2udoμy1u==Randminμy1u,maxμy1uμy2u==Randminμy2u,maxμy2uend forend forforj = 1 to 4 do (to calculate stop sign duration for a route)wju==Randminwju,maxwjuend forfori = 1 to e do (to calculate traffic light duration for a route)tiu==Rand0,maxβi,uend forforh = 1 to 5 do (to calculatez6andz7after receiving user preferences)inputσh1 and σh2calculatezhSortz6 and z7 from min to max by using Quick sort algorithmdisplay minz1−z7 values with the corresponding specific route u choicesend forend forIn 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. These constraint values are then used for calculating z2−z5. Values of σh are the inputs to the algorithm for calculating the weighted sum value z6. The values attained from z1−z5 are then used for calculating the best vector set solution z7; the decision maker can choose a route from a variety of routes by using z7 that depends on the minimum values attained for z1−z5. A real-life application of Algorithm 5 with specific entries for minimum and maximum values is given in Section 6.For instance, the setP which is an output of Algorithm 2 is a critical input to Algorithm 5. The set P consists of p^γs vectors which has the following information: (i) the geodesic γs, (ii) the length of the geodesic lγs, and (iii) the LCM circle Cγs that has γs. The stochastic parameters are calculated by using a uniform distribution.Note that running Algorithm5 for sufficient number of times and averaging the outcomes will yield the desired objective function values. The highest and lowest values of the parameters presented in constraints (15)–(18) determine the minimum and maximum values for the objective function z2. The locations and numbers of traffic signs along with the speed limits and maximum durations of traffic lights are assumed to be known.


## 6. Numerical Analysis and Simulation

In this section, we apply the methodologies from Sections4 and 5 to two different real-life routing cases, short and long distances, and analyze the results from computational experiments and discrete event simulations. The 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. The 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 short- and 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.Table 5
A summary of parameter settings.
b = 3η11=1, η21=η31=η41=0 (short-distance routing)η12=η22=η32=0, η42=3 (short-distance routing)η13=3, η23=3, η33=2, η43=1 (short-distance routing)η11=3, η21=η31=0, η41=1 (long-distance routing)η12=2, η22=η32=0, η42=1 (long-distance routing)η13=1, η23=η33=0, η43=1 (long-distance routing)wj=U∼3j,3j+3, ∀j, j = 1, 2, 3, 4I0=1δ1=16, δ2=21, δ3=11 (short-distance routing)δ1=13, δ2=15, δ3=12 (long-distance routing)tiu=U∼15,75, ∀u, u=1,…,3k = 1, 2μyku=U∼3,9, ∀u and kV = {35, 55, 75}c1u=$0.0579 per milec2u=$0.0179 per mileσh1=0.2, ∀h, h = 1, …, 5Several 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. The 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), andj-way stop signs (3j seconds, 3j + 3 seconds) where j is a type of a stop sign for all j=1,…,4 [26]. The 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 (c1u) and operational costs (c2u) 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, c1u, equal to $0.0813 per mile (i.e., $2.376/29). If a vehicle travels mostly on highways, then c1u is $0.0579 per mile (i.e., $2.376/41). To calculate the maintenance cost per mile, c2u, 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 c2u equal to $0.0179 (i.e., 242/13,476).There are three scenarios considered for each evaluation criteria of the objective functionsz1, z2, z3, and z5: best, worst, and average. The 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. The 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 perform 1000 replications of simulation of the waiting times and collect the average statistics.The best, worst, and average scenarios are different forz1; Algorithms 1–4 are used to calculate the total distance values for z1. This 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. The best case value of z1 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.Sincez3 includes the cost incurred due to waiting at traffic signals, the best- and worst-case values of z3 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, best- and worst-case values of z5 correspond to the shortest and longest duration of waiting times. Note that there are no best- and worst-case values for z4 as the number of traffic signals does not vary on a given route. Next, we analyze the computational results from the short-distance case.


### 6.1. 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. There 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.Table6 summarizes the computational results for the routes A, B, and C concerning criteria z1 to z6. The 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)
Whenz1 is considered, Route B is a good choice in all scenarios: best, worst, and average. The 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)
Whenz2 is considered, Route C is a good choice in all scenarios. Though 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. The waiting time information is not available from Google Maps. Note that the simulated average values obtained for z2 has a significant impact on total travel time.(3)
Route B is the most cost effective whenz3 is considered. The 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 z4, Route A is the best as it has a minimum number of traffic signs.(4)
Whenz5 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. The 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.(5)
Route A is a good choice with a correspondingz6 value of 0.18 when equal weights are assumed for all the criteria (i.e., wi=0.2 for all zi where i = 1, …, 5). The priority vector solution, z7, for Routes A, B, and C is the set BV = (10.84, 4.1, 0.868, 18, 15.93), (7.21, 5.65, 0.715, 29, 18.01), and (7.366, 3.3, 0.0.729, 22, 15.93). This set displays all the solution vectors and assist the decision maker(s) to examine the best options for all routes.Table 6
Multiple criteria evaluation of short-distance routes.
Evaluation criteriaScenarioRoute A (highway)Route B (local)Route C (local)z1 (miles) total travel distanceBest10.847.217.366Worst11.387.227.367Average (simulation)11.117.217.37Google Maps10.907.207.60z2 (minutes) total waiting at traffic signalsBest4.15.653.3Worst20.327.4515.4Average (simulation)12.0616.969.23Google MapsN/AN/AN/Az3 ($) total travel costBest0.8680.7150.729Worst0.9270.7380.741Average (simulation)0.9010.7270.735Google MapsN/AN/AN/Az4 (#) number of traffic signalsStop signs139Yield signs152Traffic lights162111Total182922Google MapsN/AN/AN/Az5 (minutes) total travel timeBest15.9318.0115.93Worst31.7139.8228.03Average (simulation)24.3228.5321.97Google Maps171817z6 (minutes) weighted timeBest0.180.190.19Worst0.250.250.23Overall, the above computational experiments for short-distance routes highlight the complexity of decision-making in evaluating routes when multiple criteria are considered. With criteriaz1 to z6, a commuter can pick and choose from several travel options based on individual preference. The 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.


### 6.2. Long-Distance Routing: A Case Study

We choose an origin (O) and destination (D) located in the states of Connecticut (CT) and California (CA), respectively, and restrict our attention to three different long-distance routes that connect O and D: (A) via I-80 West; (B) via I-40 West; (C) via I-70 West. Note that all the above pathways are highway routes as they connect the east and west coasts of the US.Table7 summarizes the computational results for the routes A, B, and C with respect to criteria z1 to z6. The best values are given in bold. We compare the results with Google Maps when the relevant information is available. The following are the observations extracted from Table 7:(1)
Route A is a good choice in the best-case scenario (using innermost lanes) whenz1 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. Though 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). This 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,z2, is considered. Google Maps does not provide this information; therefore, no comparisons are made. When z3 is considered, the difference between the simulation average values for all routes is minimal. However, in the long run, a small difference can make a big impact. For instance, Route B generates a savings of $0.44 per trip (i.e., 295.96 − 295.52) when compared to Route A. For a fleet of 1000 trucks and 50 trips per year, the total savings is equal to 0.44 ∗ 1000 ∗ 50 = $22,000. A similar comparison of Route B and Route C generates a savings of $208,500. With respect to z4, both Routes A and C are the good choices as they have the minimum number of traffic signals.(3)
Route A is the good choice in the best-case scenario when the total travel time,z5, is considered. However, in both worst and average cases, Route B is a good choice. Furthermore, Google Maps suggest Route A or C, while simulation averages indicate Route B is the best when the total travel time is considered. Note that the results from Google Maps are higher than our time estimates due to the difference in the speed limits that is used to calculate the travel times.(4)
When equal weights are assumed for all the criteria (i.e.,wi=0.2 for all zi where i = 1, …, 5), any route is a good choice since all of them have a corresponding z6 value of 0.16. Following is the set of priority vector solutions, z7, for Routes A, B, and C: BV = (2979.8, 3.45, 293.489, 17, 3254.23), (2986.54, 3.95, 294.18, 19, 3261.99), and (3027.41, 3.25, 298.17, 17, 3305.88). This set displays all the solution vectors and assist the decision maker(s) to examine the best options for all routes.Table 7
Multiple criteria evaluation of long-distance routes.
Evaluation criteriaScenarioRoute A (I-80)Route B (I-40)Route C (I-70)z1 (miles) total travel distanceBest2978.882986.543027.41Worst3029.333013.483058.40Average (simulation)3005.083000.033043.23Google Maps3008.003003.003043.90z2 (minutes) total waiting at it traffic signalsBest3.453.953.25Worst16.8519.3515.75Average (simulation)10.2811.859.52Google MapsN/AN/AN/Az3 ($) total travel costBest293.49294.18298.17Worst298.38296.82301.24Average (simulation)295.96295.52299.69Google MapsN/AN/AN/Az4 (#) number of traffic signalsStop signs432Yield signs332Traffic lights131512Total171917Google MapsN/AN/AN/Az5 (minutes) total travel timeBest2556.922562.972596.49Worst2609.452598.032625.4Average (simulation)2585.402582.842617.80Google Maps264027002640z6 (minutes) weighted timeBest0.160.160.16Worst0.160.160.16In 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.Though 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 Figure4. 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. Though the savings may seem minimal, when the scope is increased from one trip to hundreds of trips per year and from one truck to the entire fleet of trucks, the magnitude of the savings can be significant.Figure 4
Distance comparison—simulation vs Google Maps. (a) Short-distance routing and (b) long-distance routing.
(a)(b)


## 6.1. 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. There 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.Table6 summarizes the computational results for the routes A, B, and C concerning criteria z1 to z6. The 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)
Whenz1 is considered, Route B is a good choice in all scenarios: best, worst, and average. The 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)
Whenz2 is considered, Route C is a good choice in all scenarios. Though 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. The waiting time information is not available from Google Maps. Note that the simulated average values obtained for z2 has a significant impact on total travel time.(3)
Route B is the most cost effective whenz3 is considered. The 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 z4, Route A is the best as it has a minimum number of traffic signs.(4)
Whenz5 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. The 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.(5)
Route A is a good choice with a correspondingz6 value of 0.18 when equal weights are assumed for all the criteria (i.e., wi=0.2 for all zi where i = 1, …, 5). The priority vector solution, z7, for Routes A, B, and C is the set BV = (10.84, 4.1, 0.868, 18, 15.93), (7.21, 5.65, 0.715, 29, 18.01), and (7.366, 3.3, 0.0.729, 22, 15.93). This set displays all the solution vectors and assist the decision maker(s) to examine the best options for all routes.Table 6
Multiple criteria evaluation of short-distance routes.
Evaluation criteriaScenarioRoute A (highway)Route B (local)Route C (local)z1 (miles) total travel distanceBest10.847.217.366Worst11.387.227.367Average (simulation)11.117.217.37Google Maps10.907.207.60z2 (minutes) total waiting at traffic signalsBest4.15.653.3Worst20.327.4515.4Average (simulation)12.0616.969.23Google MapsN/AN/AN/Az3 ($) total travel costBest0.8680.7150.729Worst0.9270.7380.741Average (simulation)0.9010.7270.735Google MapsN/AN/AN/Az4 (#) number of traffic signalsStop signs139Yield signs152Traffic lights162111Total182922Google MapsN/AN/AN/Az5 (minutes) total travel timeBest15.9318.0115.93Worst31.7139.8228.03Average (simulation)24.3228.5321.97Google Maps171817z6 (minutes) weighted timeBest0.180.190.19Worst0.250.250.23Overall, the above computational experiments for short-distance routes highlight the complexity of decision-making in evaluating routes when multiple criteria are considered. With criteriaz1 to z6, a commuter can pick and choose from several travel options based on individual preference. The 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.


## 6.2. Long-Distance Routing: A Case Study

We choose an origin (O) and destination (D) located in the states of Connecticut (CT) and California (CA), respectively, and restrict our attention to three different long-distance routes that connect O and D: (A) via I-80 West; (B) via I-40 West; (C) via I-70 West. Note that all the above pathways are highway routes as they connect the east and west coasts of the US.Table7 summarizes the computational results for the routes A, B, and C with respect to criteria z1 to z6. The best values are given in bold. We compare the results with Google Maps when the relevant information is available. The following are the observations extracted from Table 7:(1)
Route A is a good choice in the best-case scenario (using innermost lanes) whenz1 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. Though 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). This 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,z2, is considered. Google Maps does not provide this information; therefore, no comparisons are made. When z3 is considered, the difference between the simulation average values for all routes is minimal. However, in the long run, a small difference can make a big impact. For instance, Route B generates a savings of $0.44 per trip (i.e., 295.96 − 295.52) when compared to Route A. For a fleet of 1000 trucks and 50 trips per year, the total savings is equal to 0.44 ∗ 1000 ∗ 50 = $22,000. A similar comparison of Route B and Route C generates a savings of $208,500. With respect to z4, both Routes A and C are the good choices as they have the minimum number of traffic signals.(3)
Route A is the good choice in the best-case scenario when the total travel time,z5, is considered. However, in both worst and average cases, Route B is a good choice. Furthermore, Google Maps suggest Route A or C, while simulation averages indicate Route B is the best when the total travel time is considered. Note that the results from Google Maps are higher than our time estimates due to the difference in the speed limits that is used to calculate the travel times.(4)
When equal weights are assumed for all the criteria (i.e.,wi=0.2 for all zi where i = 1, …, 5), any route is a good choice since all of them have a corresponding z6 value of 0.16. Following is the set of priority vector solutions, z7, for Routes A, B, and C: BV = (2979.8, 3.45, 293.489, 17, 3254.23), (2986.54, 3.95, 294.18, 19, 3261.99), and (3027.41, 3.25, 298.17, 17, 3305.88). This set displays all the solution vectors and assist the decision maker(s) to examine the best options for all routes.Table 7
Multiple criteria evaluation of long-distance routes.
Evaluation criteriaScenarioRoute A (I-80)Route B (I-40)Route C (I-70)z1 (miles) total travel distanceBest2978.882986.543027.41Worst3029.333013.483058.40Average (simulation)3005.083000.033043.23Google Maps3008.003003.003043.90z2 (minutes) total waiting at it traffic signalsBest3.453.953.25Worst16.8519.3515.75Average (simulation)10.2811.859.52Google MapsN/AN/AN/Az3 ($) total travel costBest293.49294.18298.17Worst298.38296.82301.24Average (simulation)295.96295.52299.69Google MapsN/AN/AN/Az4 (#) number of traffic signalsStop signs432Yield signs332Traffic lights131512Total171917Google MapsN/AN/AN/Az5 (minutes) total travel timeBest2556.922562.972596.49Worst2609.452598.032625.4Average (simulation)2585.402582.842617.80Google Maps264027002640z6 (minutes) weighted timeBest0.160.160.16Worst0.160.160.16In 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.Though 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 Figure4. 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. Though the savings may seem minimal, when the scope is increased from one trip to hundreds of trips per year and from one truck to the entire fleet of trucks, the magnitude of the savings can be significant.Figure 4
Distance comparison—simulation vs Google Maps. (a) Short-distance routing and (b) long-distance routing.
(a)(b)


## 7. 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. The 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. Though 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. The 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. The 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.The 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).



---

*Source: 2098495-2020-09-29.xml*