# Lagrangian Relaxation for the Multiple Constrained Robust Shortest Path Problem

Authors: Yiyong Pan
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The study focuses on a multiple constrained reliable path problem in which travel time reliability and resource constraints are collectively considered. Nonlinear optimization model is developed to model constrained robust shortest path problem. The dual nature of the proposed problem is deduced based on the Lagrangian duality theory. An efficient algorithm based on Lagrangian dual relaxation is designed to solve constrained robust shortest path problem. An extension problem that considers multiple constraints is discussed. Numerical studies indicate that the proposed algorithm is efficient in terms of obtaining the close-to-optimal solutions within reasonable computational times.


---

## Body

## 1. Introduction

### 1.1. Motivation

The constrained shortest path problem corresponds to the extension of the shortest path problem, and this is intended to obtain a shortest path that meets several constraints [1, 2]. The deficiency of classical constrained shortest path problems is that they assume the weights of each link as deterministic and not stochastic [3]. However, the weights of each link, such as link travel time, can be stochastic in reality [4–8]. Therefore, constrained shortest path problems should be extended from a deterministic network to a stochastic network.In a stochastic network, the constrained shortest path problem defines the expected value of path travel time as the objective function of path that does not consider travel time reliability [9, 10]. However, several studies indicated that travelers care about saving travel time while improving travel time reliability in the stochastic network [11–15]. The unconstrained robust shortest path problem is extensively investigated [16–25].However, to the best of the author’s knowledge, there is a lack of studies on the constrained robust shortest path problem. The author of this study made preliminary research results related to the constrained robust shortest path problem in the early stage, which only focused on the reliable shortest path problem with only one single constraint [26]. The study extends the problem to the general situation: multiple constrained robust shortest path problem. To fill the gap, the study focuses on the multiple constrained robust shortest path problem. A constrained robust shortest path problem aids in revealing the route choice behavior by considering reliability under a few resource constraints in a traffic network. However, the computational complexity of the problem significantly increases due to the nonlinear and nonadditive nature of the objective function and a few complex constraints.


### 1.2. Literature Review

The constrained shortest path problem is examined extensively. The constrained shortest path problem is an NP-hard problem and cannot be directly solved by conventional labeling algorithms [27, 28]. Classical methods for dealing with the constrained shortest path problem include the path ranking method, dynamic programming-based method, Lagrangian relaxation-based method, and heuristic method.( 1 ) Path ranking method. Handler and Zang [3] first adopted a kth shortest path algorithm to solve the constrained shortest path problem. The basic idea involves determining k shortest paths until a first feasible path that satisfies constraints is obtained. The approach is easy and exact with respect to identifying the optimal path in a small network. The disadvantage of the method is that the computational cost and storage of algorithm increase with k, and the defect causes the algorithm to be impractical in a large network.( 2 ) Dynamic programming-based method [29–35]. Modified label algorithms based on dynamic programming are used to solve the constrained shortest path problem with respect to a one side constraint or multiple side constraints. The advantage of the method is that conventional label algorithms are directly changed to solve the constrained shortest path problem, and the algorithm performs very well in solving reasonably sized networks. However, the algorithm performs very poorly in large networks due to the “curse of dimensionality” of dynamic programming.( 3 ) Lagrangian relaxation-based method [3, 9, 36–40]. The main idea of the method involves dualizing the side constraints into the objection function. The original problem is translated into a Lagrangian dual problem based on Lagrangian duality theory. The corresponding optimal value of Lagrangian dual problem is the lower bound of the primal problem. The advantage of Lagrangian relaxation-based method is that it is effective in terms of a large sized network. The disadvantage of the method is that it obtains near-optimal solutions that occasionally do not satisfy accuracy requirements and the defect is eliminated by combining this with unconstrained shortest path algorithms. Handler and Zang [3] first developed a Lagrangian relaxation-based method to solve the constrained shortest path problem and obtained near-optimal solutions. Carlyle et al. [38] presented a Lagrangian relaxation-based method combined with enumerating near-shortest paths to solve the constrained shortest path problem. Muhandiramge and Boland [39] developed an iterative algorithm in which nodes and edges were eliminated to solve the Lagrangian dual problems. Wang et al. [9] examined a stochastic constrained shortest path problem and proposed a mixed algorithm based on Lagrangian relaxation to solve it although they only focused on obtaining the least expected travel time path and did not consider travel time reliability.( 4 ) Heuristic method [34, 40–45]. In order to improve the speed of calculation in large networks, various heuristic methods are proposed. The advantage of the method is that heuristic techniques embedded in an algorithm improve the computational efficiency of the algorithm specifically in large networks. However, in the heuristic method, it is difficult to obtain an exact solution of the constrained shortest path problem, and it is occasionally difficult to even obtain a close-to-optimal solution. The heuristic technique is only used as an auxiliary means to improve the computational efficiency of an algorithm.The reliable shortest path problem is examined by several researchers. It is formulated in various forms based on the objection function of path including the(1) most reliable path problem (Frank, 1969; Nie and Wu, 2009a), (2)  α-reliable path problem (Chen and Ji, 2005; Chen et al., 2012; [19]), and (3) mean-standard deviation shortest path problem [18, 20, 21, 23, 46, 47]. The first two problems are complementary to each other. The α-reliable path problem is translated into the mean-standard deviation shortest path problem when travel times are normally distributed. The third problem defines the objective function of path as the sum of mean and standard deviation of stochastic travel time, and they are only constructed based on the stochastic parameter that efficiently reduces the computation cost and memory usage. The study focuses on the third problem with resource constraints and extends studies where the third problem is discussed.Sen et al. [46] considered the reliable path problem by defining the objection function as the sum of a mean and a variance of path. The solution to the problem is relatively simple due to the additive nature of the objective function although its disadvantage is that the mean and variance are not in an order of magnitude. Xing and Zhou [18] considered the mean-standard deviation path problem that directly reflects the reliability of the path although the solution to the problem is complicated due to the nonlinear and nonadditive nature of the objective function. A Lagrangian dual relaxation algorithm is constructed to solve the problem. Zeng et al. [47] focused on the reliable path problem in a stochastic network with spatial correlation. The Cholesky decomposition is used to resolve the correlation of link travel time, and a similar Lagrangian dual relaxation algorithm is developed to solve it. Khani and Boyles [20] discussed the relationship between the mean-standard deviation path problem and the mean-variance path problem and developed an iterative algorithm to solve the mean-standard deviation path problem. Wu [21] extended the mean-standard deviation path problem from a stochastic and static network to a stochastic and time-dependent network and introduced stochastic dominance theory to solve the problem. Zhang et al. [23] reconstructed the mean-standard deviation path problem to a convex problem. They analyzed duality gaps in different Lagrangian relaxation approaches and selected a novel Lagrangian relaxation algorithm to solve it.


### 1.3. Contribution

The study first focuses on a constrained robust shortest path problem, and the specific contributions are as follows:( 1 ) The study first proposes a multiple constrained robust shortest path problem to obtain the mean-standard deviation shortest path subject to a few resource constraints in stochastic network. A reliable objective function and a few resource constraints are simultaneously considered. The constrained robust shortest path problem is formulated as a mixed 0-1 convex nonlinear programming model.( 2 ) The dual problem of constrained robust shortest path problem is derived based on Lagrangian duality theory, and the dual problem is divided into two independent subproblems that are easy to solve.( 3 ) An efficient algorithm is constructed to solve a constrained robust shortest path problem. A gradient descent algorithm is developed to solve the dual problem by gradually reducing the gap between the upper and lower bounds. A kth shortest path algorithm is used to obtain the feasible path.( 4 ) An extension problem that considers multiple constraints is further discussed, and the corresponding dual problem is also derived.( 5 ) A numerical test is developed on a small network, medium-scale transportation network, and large-scale network. The numerical results are analyzed to verify the effectiveness and stability of the proposed algorithm.The study is organized as follows. In Section2, the constrained robust shortest path problem is discussed and formulated as a mixed 0-1 convex nonlinear programming model. In Section 3, the dual problem is derived. In Section 4, an efficient algorithm is constructed. In Section 5, an extension problem that considers multiple constraints is discussed. Numerical experiments are elucidated in Section 6. Finally, Section 7 details a few conclusions and discusses future research directions.


## 1.1. Motivation

The constrained shortest path problem corresponds to the extension of the shortest path problem, and this is intended to obtain a shortest path that meets several constraints [1, 2]. The deficiency of classical constrained shortest path problems is that they assume the weights of each link as deterministic and not stochastic [3]. However, the weights of each link, such as link travel time, can be stochastic in reality [4–8]. Therefore, constrained shortest path problems should be extended from a deterministic network to a stochastic network.In a stochastic network, the constrained shortest path problem defines the expected value of path travel time as the objective function of path that does not consider travel time reliability [9, 10]. However, several studies indicated that travelers care about saving travel time while improving travel time reliability in the stochastic network [11–15]. The unconstrained robust shortest path problem is extensively investigated [16–25].However, to the best of the author’s knowledge, there is a lack of studies on the constrained robust shortest path problem. The author of this study made preliminary research results related to the constrained robust shortest path problem in the early stage, which only focused on the reliable shortest path problem with only one single constraint [26]. The study extends the problem to the general situation: multiple constrained robust shortest path problem. To fill the gap, the study focuses on the multiple constrained robust shortest path problem. A constrained robust shortest path problem aids in revealing the route choice behavior by considering reliability under a few resource constraints in a traffic network. However, the computational complexity of the problem significantly increases due to the nonlinear and nonadditive nature of the objective function and a few complex constraints.


## 1.2. Literature Review

The constrained shortest path problem is examined extensively. The constrained shortest path problem is an NP-hard problem and cannot be directly solved by conventional labeling algorithms [27, 28]. Classical methods for dealing with the constrained shortest path problem include the path ranking method, dynamic programming-based method, Lagrangian relaxation-based method, and heuristic method.( 1 ) Path ranking method. Handler and Zang [3] first adopted a kth shortest path algorithm to solve the constrained shortest path problem. The basic idea involves determining k shortest paths until a first feasible path that satisfies constraints is obtained. The approach is easy and exact with respect to identifying the optimal path in a small network. The disadvantage of the method is that the computational cost and storage of algorithm increase with k, and the defect causes the algorithm to be impractical in a large network.( 2 ) Dynamic programming-based method [29–35]. Modified label algorithms based on dynamic programming are used to solve the constrained shortest path problem with respect to a one side constraint or multiple side constraints. The advantage of the method is that conventional label algorithms are directly changed to solve the constrained shortest path problem, and the algorithm performs very well in solving reasonably sized networks. However, the algorithm performs very poorly in large networks due to the “curse of dimensionality” of dynamic programming.( 3 ) Lagrangian relaxation-based method [3, 9, 36–40]. The main idea of the method involves dualizing the side constraints into the objection function. The original problem is translated into a Lagrangian dual problem based on Lagrangian duality theory. The corresponding optimal value of Lagrangian dual problem is the lower bound of the primal problem. The advantage of Lagrangian relaxation-based method is that it is effective in terms of a large sized network. The disadvantage of the method is that it obtains near-optimal solutions that occasionally do not satisfy accuracy requirements and the defect is eliminated by combining this with unconstrained shortest path algorithms. Handler and Zang [3] first developed a Lagrangian relaxation-based method to solve the constrained shortest path problem and obtained near-optimal solutions. Carlyle et al. [38] presented a Lagrangian relaxation-based method combined with enumerating near-shortest paths to solve the constrained shortest path problem. Muhandiramge and Boland [39] developed an iterative algorithm in which nodes and edges were eliminated to solve the Lagrangian dual problems. Wang et al. [9] examined a stochastic constrained shortest path problem and proposed a mixed algorithm based on Lagrangian relaxation to solve it although they only focused on obtaining the least expected travel time path and did not consider travel time reliability.( 4 ) Heuristic method [34, 40–45]. In order to improve the speed of calculation in large networks, various heuristic methods are proposed. The advantage of the method is that heuristic techniques embedded in an algorithm improve the computational efficiency of the algorithm specifically in large networks. However, in the heuristic method, it is difficult to obtain an exact solution of the constrained shortest path problem, and it is occasionally difficult to even obtain a close-to-optimal solution. The heuristic technique is only used as an auxiliary means to improve the computational efficiency of an algorithm.The reliable shortest path problem is examined by several researchers. It is formulated in various forms based on the objection function of path including the(1) most reliable path problem (Frank, 1969; Nie and Wu, 2009a), (2)  α-reliable path problem (Chen and Ji, 2005; Chen et al., 2012; [19]), and (3) mean-standard deviation shortest path problem [18, 20, 21, 23, 46, 47]. The first two problems are complementary to each other. The α-reliable path problem is translated into the mean-standard deviation shortest path problem when travel times are normally distributed. The third problem defines the objective function of path as the sum of mean and standard deviation of stochastic travel time, and they are only constructed based on the stochastic parameter that efficiently reduces the computation cost and memory usage. The study focuses on the third problem with resource constraints and extends studies where the third problem is discussed.Sen et al. [46] considered the reliable path problem by defining the objection function as the sum of a mean and a variance of path. The solution to the problem is relatively simple due to the additive nature of the objective function although its disadvantage is that the mean and variance are not in an order of magnitude. Xing and Zhou [18] considered the mean-standard deviation path problem that directly reflects the reliability of the path although the solution to the problem is complicated due to the nonlinear and nonadditive nature of the objective function. A Lagrangian dual relaxation algorithm is constructed to solve the problem. Zeng et al. [47] focused on the reliable path problem in a stochastic network with spatial correlation. The Cholesky decomposition is used to resolve the correlation of link travel time, and a similar Lagrangian dual relaxation algorithm is developed to solve it. Khani and Boyles [20] discussed the relationship between the mean-standard deviation path problem and the mean-variance path problem and developed an iterative algorithm to solve the mean-standard deviation path problem. Wu [21] extended the mean-standard deviation path problem from a stochastic and static network to a stochastic and time-dependent network and introduced stochastic dominance theory to solve the problem. Zhang et al. [23] reconstructed the mean-standard deviation path problem to a convex problem. They analyzed duality gaps in different Lagrangian relaxation approaches and selected a novel Lagrangian relaxation algorithm to solve it.


## 1.3. Contribution

The study first focuses on a constrained robust shortest path problem, and the specific contributions are as follows:( 1 ) The study first proposes a multiple constrained robust shortest path problem to obtain the mean-standard deviation shortest path subject to a few resource constraints in stochastic network. A reliable objective function and a few resource constraints are simultaneously considered. The constrained robust shortest path problem is formulated as a mixed 0-1 convex nonlinear programming model.( 2 ) The dual problem of constrained robust shortest path problem is derived based on Lagrangian duality theory, and the dual problem is divided into two independent subproblems that are easy to solve.( 3 ) An efficient algorithm is constructed to solve a constrained robust shortest path problem. A gradient descent algorithm is developed to solve the dual problem by gradually reducing the gap between the upper and lower bounds. A kth shortest path algorithm is used to obtain the feasible path.( 4 ) An extension problem that considers multiple constraints is further discussed, and the corresponding dual problem is also derived.( 5 ) A numerical test is developed on a small network, medium-scale transportation network, and large-scale network. The numerical results are analyzed to verify the effectiveness and stability of the proposed algorithm.The study is organized as follows. In Section2, the constrained robust shortest path problem is discussed and formulated as a mixed 0-1 convex nonlinear programming model. In Section 3, the dual problem is derived. In Section 4, an efficient algorithm is constructed. In Section 5, an extension problem that considers multiple constraints is discussed. Numerical experiments are elucidated in Section 6. Finally, Section 7 details a few conclusions and discusses future research directions.


## 2. Problem Statement

We consider a directed graphG=(N,A) that consists of a set of nodes N=1,2,⋯,n(N=n) and a set of links A=(i,j)∣i,j∈N,i≠j(A=m). The travel time on each link (i,j) is a random variable tij, and its mean and standard deviation are cij and σij, respectively. Additionally, each link possesses a resource weight wij. A link travel time distribution is independent.The travel time along pathx is also a random variable as follows:(1)tx=∑i,j∈Atijxij,where(2)xij=1,iflinki,jbelongstothesolutionpathx0,otherwise.The path is denoted as a vector x=xij,∀i,j∈A. The mean and standard deviation of path travel time are as follows:(3)Etx=∑i,j∈AEtijxij=∑i,j∈Acijxij,(4)σtx=∑i,j∈Aσij2xij.The objective function of path travel time is given as follows:(5)fx=∑i,j∈Acijxij+λ∑i,j∈Aσij2xij,where λ denotes the reliability coefficient. The resource consumption along path x is also a variable as follows:(6)wx=∑i,j∈Awijxij.The resource consumption along pathwx must be bound by constraint conditions (e.g., fuel consumption or power consumption). Subsequently, the following expression is obtained:(7)∑i,j∈Awijxij≤W.The parameter W denotes the upper limit value allowed for resource consumption along path.Given the definition, the constrained robust shortest path problem is described as the following mathematical model:(8)min∑i,j∈Acijxij+λ∑i,j∈Aσij2xijs.t.:(9)∑i,j∈Axij-∑j,i∈Axji=1,i=O0,i∈N-O,D-1,i=D(10)∑i,j∈Awijxij≤W(11)xij∈0,1,∀i,j∈A.The constrained robust shortest path problem is more complex than the constrained shortest path problem due to the addition of the standard deviation in the objective function (8) that completely corresponds to a nonadditive function.


## 3. The Dual Problem

### 3.1. Model Reformulation

It is difficult to solve the constrained robust shortest path problem due to the nonlinear objective function (8) and constraint (9). The objective function (8) is a nonlinear and nonadditive function due to the standard deviation part. Our goal involves translating the nonadditive objective function into an additive objective function. We introduce artificial variable y to eliminate the nonadditivity of the standard deviation in the objective function.The objective function (8) is equivalent to the following functions:(12)∑i,j∈Acijxij+λy,∑i,j∈Aσij2xij=y.Lemma 1 (see [18]).
0 ≤ y ≤ y ′, where y′ denotes the variance of the least expected travel time path in network with respect to the link weight tij.Subsequently, the problem ((8)-(10)) is reformulated as the following constrained shortest path problem.(13)min∑i,j∈Acijxij+λys.t.:(14)∑i,j∈Aσij2xij≤y(15)∑i,j∈Awijxij≤W(16)∑i,j∈Axij-∑j,i∈Axji=1,i=O0,i∈N-O,D-1,i=D(17)0≤y≤y′xij∈0,1,∀i,j∈A


### 3.2. Lagrangian Dual Problem

Based on the Lagrangian duality theory, two hard constraints (14) and (15) are dualized into the objective function with two Lagrangian multipliers. Therefore, with respect to any Lagrangian multiplier vector μ=(μ1,μ2)≥0 and x∈X-,y∈Y-, where X-=x∣∑i,j∈Axij-∑j,i∈Axji=b,xij∈0,1,∀(i,j)∈A, (18)b=1,i=O0,i∈N-O,D-1,i=D.,whereY-=y0≤y≤y′.The Lagrangian function is defined as follows:(19)Lμ,x,y=∑i,j∈Acijxij+λy+μ1∑i,j∈Aσij2xij-y+μ2∑i,j∈Awijxij-W,where (20)Lμ,x=∑i,j∈Acij+μ1σij2+μ2wijxij,(21)Lμ,y=λy-μ1y-μ2W,(22)Lμ=minx∈X-,y∈Y-⁡Lμ,x,y=minx∈X-⁡Lμ,x+miny∈Y-⁡Lμ,y.The first subproblem is expressed as follows:(23)minx∈X-⁡Lμ,x.It corresponds to the classical shortest path problem with a new weight cij+μ1σij2+μ2wij. Dijkstra’s algorithm (or Dijkstra’s Shortest Path First algorithm, SPF algorithm) is an algorithm for finding the classical shortest paths between nodes in a graph [28].The second subproblem is expressed as follows:(24)miny∈Y-⁡Lμ,y.This is a nonlinear concave optimization problem, and it is easily solved through minimizing the concave function value at two interval endpoints: 0 and y′ [32]; thus the optimal solution is min⁡-μ2W,λy′-μ1y′-μ2W.Based on the Lagrangian duality theory [48], the following expression is obtained:(25)Lμ≤fx∗.The Lagrangian dual problem of the constrained robust shortest path problem ((8)-(10)) is as follows:(26)Lμ∗=maxμ≥0⁡Lμ≤fx∗.


## 3.1. Model Reformulation

It is difficult to solve the constrained robust shortest path problem due to the nonlinear objective function (8) and constraint (9). The objective function (8) is a nonlinear and nonadditive function due to the standard deviation part. Our goal involves translating the nonadditive objective function into an additive objective function. We introduce artificial variable y to eliminate the nonadditivity of the standard deviation in the objective function.The objective function (8) is equivalent to the following functions:(12)∑i,j∈Acijxij+λy,∑i,j∈Aσij2xij=y.Lemma 1 (see [18]).
0 ≤ y ≤ y ′, where y′ denotes the variance of the least expected travel time path in network with respect to the link weight tij.Subsequently, the problem ((8)-(10)) is reformulated as the following constrained shortest path problem.(13)min∑i,j∈Acijxij+λys.t.:(14)∑i,j∈Aσij2xij≤y(15)∑i,j∈Awijxij≤W(16)∑i,j∈Axij-∑j,i∈Axji=1,i=O0,i∈N-O,D-1,i=D(17)0≤y≤y′xij∈0,1,∀i,j∈A


## 3.2. Lagrangian Dual Problem

Based on the Lagrangian duality theory, two hard constraints (14) and (15) are dualized into the objective function with two Lagrangian multipliers. Therefore, with respect to any Lagrangian multiplier vector μ=(μ1,μ2)≥0 and x∈X-,y∈Y-, where X-=x∣∑i,j∈Axij-∑j,i∈Axji=b,xij∈0,1,∀(i,j)∈A, (18)b=1,i=O0,i∈N-O,D-1,i=D.,whereY-=y0≤y≤y′.The Lagrangian function is defined as follows:(19)Lμ,x,y=∑i,j∈Acijxij+λy+μ1∑i,j∈Aσij2xij-y+μ2∑i,j∈Awijxij-W,where (20)Lμ,x=∑i,j∈Acij+μ1σij2+μ2wijxij,(21)Lμ,y=λy-μ1y-μ2W,(22)Lμ=minx∈X-,y∈Y-⁡Lμ,x,y=minx∈X-⁡Lμ,x+miny∈Y-⁡Lμ,y.The first subproblem is expressed as follows:(23)minx∈X-⁡Lμ,x.It corresponds to the classical shortest path problem with a new weight cij+μ1σij2+μ2wij. Dijkstra’s algorithm (or Dijkstra’s Shortest Path First algorithm, SPF algorithm) is an algorithm for finding the classical shortest paths between nodes in a graph [28].The second subproblem is expressed as follows:(24)miny∈Y-⁡Lμ,y.This is a nonlinear concave optimization problem, and it is easily solved through minimizing the concave function value at two interval endpoints: 0 and y′ [32]; thus the optimal solution is min⁡-μ2W,λy′-μ1y′-μ2W.Based on the Lagrangian duality theory [48], the following expression is obtained:(25)Lμ≤fx∗.The Lagrangian dual problem of the constrained robust shortest path problem ((8)-(10)) is as follows:(26)Lμ∗=maxμ≥0⁡Lμ≤fx∗.


## 4. Solution Algorithm

An efficient algorithm is designed to solve constrained robust shortest path problem based on the algorithm suggested by Fisher [48]. A gradient descent iterative algorithm is proposed to solve the dual problem. Specifically, in each iteration, the first subproblem (23) and second subproblem (24) are solved, and the sum of two corresponding optimal values is set as the lower bound of the optimal value of the original problem ((8)-(10)). Furthermore, the solution of the first subproblem (23) provides a potentially feasible solution for the original problem ((8)-(10)), and the corresponding objective value of the original problem ((8)-(10)) is set as the upper bound of the optimal value of the original problem ((8)-(10)). However, the solution of the first subproblem (23) potentially does not satisfy the constraints and is not a feasible solution to the original problem ((8)-(10)). Therefore, the k shortest path algorithm is introduced to obtain a feasible solution of the original problem ((8)-(10)). The k shortest path algorithm is a generalization of the shortest path algorithm. The algorithm finds not only the shortest path, but also k − 1 other paths in nondecreasing order of cost. k is the number of shortest paths to find. Since 1957 there have been many papers published on the k shortest path routing algorithm problem. Handler and Zang [3] first adopted the k shortest path algorithm to solve the constrained shortest path problem. The basic idea involves determining k shortest paths until a first feasible path that satisfies constraints is obtained. Since the k shortest path algorithm itself can solve the constrained path problem, the proposed algorithm adopts this algorithm to accelerate and guarantee the convergence of the proposed algorithm. In general, the upper bound and lower bound of the proposed algorithm change during the iterative process. In extreme cases, the upper bound of the proposed algorithm remains unchanged during the iterative process if the upper bound of the original problem ((8)-(10)) (that is, the optimal solution) is obtained by the k shortest path algorithm, but the lower bound approximates the upper bound, which is the optimal value.In the iterative process, the gap between the upper and lower bounds gradually decreases. Additionally, when the gap corresponds to the minimum gap or zero, the corresponding upper bound is a near-optimal solution or an accurate solution of the original problem ((8)-(10)). The detailed procedure of the algorithm is described in Algorithm 1.Algorithm 1:Subgradient method for updating Lagrangian multipliers.
1. Initialization Set initial iterationk=1, initial Lagrangian multipliers μk=(μ1k,μ2k)∈R2, initial lower boundLk=-∞, initial upper bound Uk=+∞, maximum relative gap ε^, maximum value of iterationK. 2. Update the lower bound Solveminx∈X-⁡L(μk,x) obtain an optimal solution xk and optimal value Lxk. Solveminy∈Y-⁡L(μk,y) obtain an optimal solution yk and optimal value Lyk. Set the lower boundLk=max⁡Lk,Lxk+Lyk 3. Update the upper bound Ifxk satisfies all the side constraints (9) then set the upper bound: Uk=min⁡Uk,∑(i,j)∈Acijxijk+λ∑(i,j)∈Aσij2xijk. Else if, solve the problem ((8)-(10)) by K-shortest path algorithm to obtain a feasible Solutionx′k, and set the upper bound: Uk=min⁡Uk,∑(i,j)∈Acijxij′k+λ∑(i,j)∈Aσij2xij′k. End 4. Update Lagrangian multipliers Compute sub-gradient direction ∇Lμ1(μ)=∑(i,j)∈Aσij2xijk-yk ∇Lμ2(μ)=∑(i,j)∈Awijxijk-W Updated Lagrangian multipliers μ1k=max⁡0,μ1k+γk1∇Lμ1(μ) μ2k=max⁡0,μ2k+γk2∇Lμ2(μ) whereγk1=δk1Uk-Lk/∑(i,j)∈Aσij2xijk-yk2, γk2=δk2Uk-Lk/∑(i,j)∈Awijsxij-Ws2, 0 ≤ δ k 1 , δ k 2 ≤ 2 the value of δk1,δk2 is suggested by Fisher [48]. 5. Convergence test If relative gapεk=Uk-Lk/Uk×100≤ε^ or k≥K stop; Otherwise k=k+1, go to step 2.


## 5. Extension to the Multiple Constraints

In this section, the reliable shortest path problem with a single constraint is extended to the problem with multiple constraints. In a stochastic network, each link possesses multiple resource weights,wijs,(s=1,2⋯,S). Given the definition, the constrained robust shortest path problem is described as follows.(27)min∑i,j∈Acijxij+λ∑i,j∈Aσij2xijs.t.:(28)∑i,j∈Awijsxij≤Ws,s=1,2⋯,S(29)x∈X-The parameterWs denotes the upper value of the resource weights. The constrained robust shortest path problem ((27)-(29)) is reformulated as the following constrained shortest path problem:(30)min∑i,j∈Acijxij+λys.t.:(31)∑i,j∈Aσij2xij≤y(32)∑i,j∈Awijsxij≤Ws,s=1,2⋯,S(33)x∈X-,y∈Y-where X-=x∣∑i,j∈Axij-∑j,i∈Axji=b,xij∈0,1,∀(i,j)∈A, Y-=y∣0≤y≤y′.For any Lagrangian multiplier vectorμ=(μ1,μ2⋯,μs+1)≥0 and x∈X-,y∈Y-, the Lagrangian function is defined as follows:(34)Lμ,x,y=∑i,j∈Acijxij+λy+μ1∑i,j∈Aσij2xij-y+∑s=1Sμs+1∑i,j∈Awijsxij-Wswhere (35)Lμ,x=∑i,j∈Acij+μ1σij2+∑s=1Sμs+1wijsxij(36)Lμ,y=λy-μ1y-∑s=1Sμs+1Ws,(37)Lμ=minx∈X-,y∈Y-⁡Lμ,x,y=minx∈X-⁡Lμ,x+miny∈Y-⁡Lμ,y.The Lagrangian dual problem of the constrained robust shortest path problem ((30)-(33)) is as follows:(38)Lμ∗=maxμ≥0⁡Lμ≤fx∗.A solution algorithm suggested in Section4 is modified to solve the multiple constrained robust shortest path problem. Step 4 in the algorithm updating Lagrangian multipliers is modified as follows, and the other steps are similar to the solution algorithm suggested in Section 4.Compute the subgradient direction(39)∇Lμ1μ=∑i,j∈Aσij2xijk-yk(40)∇Lμs+1μ=∑i,j∈Awijsxij-Ws,s=1,2⋯,SUpdate Lagrangian multipliers(41)μ1k=max⁡0,μ1k+γk1∇Lμ1μ(42)μs+1k=max⁡0,μs+1k+γk2∇Lμs+1μ,s=1,2⋯,Swhere γk1=δk1Uk-Lk/∑(i,j)∈Aσij2xijk-yk2, γk2=δk2Uk-Lk/∑s=1S∑(i,j)∈Awijsxij-Ws2,0≤δk1,δk2≤2, and the values of δk1,δk2 are suggested by Fisher [48].


## 6. Numerical Experiments

This section verifies the validity and feasibility of the proposed model and algorithm. All algorithms were coded via Matlab and tested on a Windows-10 (64) workstation with two 2.00 GHz Xeon CPUs and 4G RAM. Three networks are considered to represent different experimental environments and include a small network, Sioux Falls network [49], and large-scale network [6].


### 6.1. A Small Network

A small network is constructed with 5 nodes and 6 links as shown in Figure2. Link travel time tij is a random variable with its mean cij and standard deviation σij. Additionally, each link possesses two resource weights wij1 and wij2. The parameters of each link [cij,σij2,wij1,wij2] are shown in Figure 1.Figure 1
Small network.Figure 2
Updated process of the upper bound, lower bound, and relation gap with respect to the iterations.The initial parameters are set as the OD pair:1-5,ε^=1,λ=0.5, and K=100. The cases of different constraints are computed. Table 1 shows the updated process of Lagrangian multiplier, bound, and gap based on iterations with two upper values of the resource weights: W1=10 and W2=20. Figure 3 shows the updated process of the bound and relative gaps with respect to the iteration. Table 2 shows the updated process of the bound and gap based on iterations with W1=10 and W2=10. Table 3 shows the updated process of the bound and gap based on iterations with W1=20 and W2=20.Table 1
Updated process of the Lagrangian multiplier, bound, and gap with respect to the iterations (W1=10,W2=20).
k μ 1 k μ 2 k μ 3 k U k L k U k - L k ε k 1 100.0 100.0 100.0 60.00 0.00 60.00 100.0  2 61.31 61.48 62.47 52.26 27.20 25.06 47.95  3 63.49 68.23 61.10 48.20 31.12 17.08 35.43  4 53.49 66.15 51.29 47.40 35.49 11.91 25.12  5 51.43 59.30 53.20 43.26 35.98 7.28 16.82  6 30.20 35.36 50.20 42.16 38.32 3.84 9.10  7 34.23 0.00 50.20 40.00 39.13 0.87 2.17  8 32.45 0.00 50.20 40.00 40.00 0.00 0.00Table 2
Updated process of the bound and gap with respect to the iterations (W1=10,W2=10).
k U k L k U k - L k ε k 1 60.00 0.00 60.00 100  2 58.26 29.20 29.06 49.87  3 48.93 36.82 12.11 24.74  4 47.73 45.49 2.24 4.69  5 47.73 47.73 0.00 0.00Table 3
Updated process of the bound and gap with respect to the iterations (W1=20,W2=20).
k U k L k U k - L k ε k 1 60.00 0.00 60.00 100.00  2 52.75 23.20 29.55 56.01  3 45.20 26.12 19.08 42.21  4 43.53 25.40 18.13 41.64  5 38.81 34.33 4.48 11.54  6 38.87 38.87 0.00 0.00Figure 3
Sioux Falls network.As shown in Table1, the exact optimal value is 40, obtained at the 8th iteration. The two resource weights of optimal path 1-3-4-5 are 9 and 12, and this satisfies two constraints. As shown in Figure 2, the upper bounds are gradually close to the lower bounds, and an exact optimal value is finally determined when the relation gap is reduced to 0 in the 8th iteration. As shown in Table 2, an exact optimal value is 47.73, and this is obtained at the 5th iteration. The two resource weights of optimal path 1-2-5 are 7 and 9, and this satisfies two constraints. As shown in Table 3, an exact optimal value is 38.87, and this is obtained at the 6th iteration. The two resource weights of the optimal path 1-2-4-5 are 11 and 12, and this satisfies two constraints. There are differences in the optimal value and optimal path with respect to the constraint conditions.


### 6.2. Sioux Falls Network

Sioux Falls network is a medium-scale transportation network that includes 24 zones, 24 nodes, and 76 links as shown in Figure3 [49].


#### 6.2.1. Experimental Design

The experimental design is set as follows. In the Sioux Falls network, link travel timestij((i,j)∈A) are random variables, its mean is cij((i,j)∈A), and its standard deviation is σij((i,j)∈A). In order to simulate its stochastic properties, the expected values and variances are randomly obtained from certain intervals as shown in Table 4. Each link possesses a resource weight wij that is also randomly obtained from a certain interval.Table 4
Intervals of randomly generated link parameters for different links.
c i j σ i j w i j ( i , j ) [ 1,18 ] [ 1,25 ] [ 1,30 ] ( 3,1 ),(3,4),(3,12),(9,5),(9,8),(9,10),(24,13),(24,21),(24,23)  [ 1,25 ] [ 1,25 ] [ 1,20 ] ( 5,4 ),(5,6),(5,9),(10,9),(10,11),(10,15),(10,16),(10,17)  [ 1,20 ] [ 1,28 ] [ 1,25 ] ( 15,10 ),(15,14),(15,19),(15,22),(20,18),(20,19),(20,21),(20,22)  [ 1,18 ] [ 1,20 ] [ 1,19 ] ( 4,3 ),(4,5),(4,11),(8,6),(8,7),(8,9),(8,16),(12,3),(12,11),(12,13),(16,8),(16,10),(16,17),(16,18),  [ 1,15 ] [ 1,30 ] [ 1,20 ] ( 7,8 ),(7,18),(14,11),(14,15),(14,23),(21,20),(21,22),(21,24)  [ 1,20 ] [ 1,20 ] [ 1,20 ] others


#### 6.2.2. Computational Results of Different OD Pairs

Table5 shows the computational results of different OD pairs in which the initial parameters are set as ε^=5, W=50, λ=0.2, K=100. The 10 OD pairs are randomly selected from all 24 nodes. As shown in Table 5, all relative gaps are less than the given value ε^=5%, and this implies that the algorithms are terminated within the maximum step K=100. Two relative gaps are zero, and this implies that the accurate solution is finally obtained. All path lengths are less than the constraint value W=50 and satisfy the path length constraint, and this verifies the feasibility and accuracy of the proposed algorithm.Table 5
Computational results of the upper bound, lower bound, and relation gap with respect to the different OD pairs.
OD U L ε Optimal path 1 - 3 41.32 40.45 2.10 1 - 3 - 12 - 11 - 14 - 23 - 24 - 12 - 11 - 14 - 23 - 24  1 - 24 61.94 58.87 4.95 1 - 24  1 - 24 57.31 55.6 2.98 1 - 24  1 - 24 56.43 55.5 1.64 1 - 24  1 - 24 59.24 57.91 2.24 1 - 24  1 - 24 51.66 50.02 3.17 1 - 24  1 - 24 40.60 40.60 0.00 1 - 24  1 - 24 49.13 46.89 4.55 1 - 24  1 - 24 49.08 49.08 0.00 1 - 24  1 - 24 52.14 51.23 1.74 1 - 24


#### 6.2.3. Computational Results of Different Constraints

Tables6–8 show the computational results of different constraints with the same OD pair as follows:1-24, 7-19, and 2-23, respectively. The notation NO in Tables 6 and 7 indicates that the algorithm is not successful in obtaining the optimal value when the constrained value W is less than the shortest path weight, a path that satisfies the constraint is absent, and the constrained robust shortest path problem does not include a feasible solution. In Tables 6–8, W=+∞ implies that the constrained robust shortest path problem is reduced to the reliable shortest path problem without constraints. The optimal path varies with respect to the constraint condition that is denoted by constraint value. For example, in Table 6, the optimal paths are different when the constraint values are W=50 and60. Similarly, in Table 7, the optimal paths are also different when the constraint values are W=30 and40. Finally, in Table 8, the optimal paths are also different when the constraint values are W=40,50, and55. Specifically, differences exist in the optimal paths solved by the reliable shortest path problem with and without constraint. For example, in Table 8, the optimal paths are different when the constraint values are W=55 and +∞. Conversely, with increases in the constraint value and path length, the upper bound and lower bound of the optimal path exhibit a decreased trend. For example, in Table 8, the constraint value increases from 40 to 55 and the corresponding path length increases from 36 to 54 although the upper bound decreases from 72.44 to 57.12 and the lower bound decreases from 70.39 to 54.39. The relative gap is not significantly influenced by the constraint condition.Table 6
Computational results of the upper bound, lower bound, and relation gap with respect to the different constraints (1-24).
W U L ε Optimal path 20 NO NO NO NO  50 41.32 40.45 2.10 1 - 3 - 12 - 11 - 14 - 23 - 24  60 37.56 36.12 3.83 + ∞  + ∞ 37.56 37.56 0.00 7 - 19Table 7
Computational results of the upper bound, lower bound, and relation gap with respect to the different constraints (7-19).
W U L ε Optimal path 20 NO NO NO NO  30 51.31 51.48 -0.03 7 - 18 - 16 - 17 - 19  40 33.49 33.49 0.00 7 - 18 - 20 - 19  + ∞ 33.49 33.49 0.00 7 - 18 - 20 - 19Table 8
Computational results of the upper bound, lower bound, and relation gap with respect to the different constraints (2-23).
W U L ε Optimal path 40 72.44 70.39 2.83 2 - 6 - 5 - 4 - 11 - 14 - 23  50 61.94 58.87 4.95 55  55 57.12 54.39 4.77 2 - 1 - 3 - 12 - 11 - 14 - 23  + ∞ 53.20 53.20 0.00 2 - 1 - 3 - 12 - 13 - 24 - 23


#### 6.2.4. Computational Results of Different Reliability Coefficients

Table9 shows the computational results of different reliability coefficients in which the initial parameters are set as the OD pair as follows: 14-20, ε^=5, W=50, and K=50. With respect to the first case in Table 9, λ=0 implies that the constrained robust shortest path problem is reduced to the constrained shortest path problem and differences exist between the optimal paths solved by the constrained robust shortest path problem and constrained shortest path problem. The optimal path varies with respect to the reliability coefficient that denotes travelers’ attitudes to path reliability. For example, in Table 9, the optimal paths are different when the reliability coefficient λ=0.2, 0.5, and 0.8.Table 9
Computational results of the upper bound, lower bound, and relation gap with respect to the different reliability coefficients.
λ U L ε Optimal path 0.0 31.31 31.31 0.00 14 - 15 - 19 - 20  0.2 35.12 35.01 0.31 14 - 15 - 19 - 20  0.5 37.25 37.25 0.00 14 - 15 - 22 - 20  0.8 53.45 51.16 4.28 14 - 23 - 22 - 21 - 20


### 6.3. Large-Scale Network

#### 6.3.1. Experimental Design

The grid network is a scalable network that is typically used to test traffic network equilibrium problem [6]. A grid network of level h contains N=h2 nodes and A=2h(h-1) links as shown in Figure 4: h=4. The experimental design is set as follows. In the grid network, link travel times tij((i,j)∈A) are random variables, its mean is cij((i,j)∈A), and its standard deviation is σij((i,j)∈A). In order to simulate its stochastic properties, all means and standard variances are randomly selected from the interval [0,20]. Each link possesses a resource weight wij that is also randomly selected from the interval [0,20].Figure 4
Grid network.The length constraint of pathW is randomly selected from the interval [0,40×h]. The initial parameters are set as ε^=5, W=50, K=50. The computational times of the algorithm are analyzed under different network scales. All CPU times were measured in seconds with the utility routine CPU time.


#### 6.3.2. Computational Time of the Proposed Algorithm

Table10 shows the average computational time of the proposed algorithm for different network scales. The network scale is represented by the level, number of nodes, and links of the grid network. It should be noted that the relative gap ε is the average of 10 computational results with the same OD pairs. The number of iterations k and CPU time is also the average of 10 computational results for the same OD pairs.Table 10
Computational results of different reliability coefficients.
Network scale ε k CPU time (sec.) k N A 5 25 40 0.00 7 4  10 100 180 0.78 11 19  20 400 760 2.41 23 52  30 900 1740 2.23 21 96  50 2500 4900 4.28 41 291  70 4900 9600 4.02 43 485  90 8100 16020 3.19 59 805Evidently, as shown in Table10, the relative gaps are all less than the initial parameters ε^=5. The CPU time significantly increases with respect to the network scale. For example, the CPU time increases from 4 S to 805 S when the level of grid network increases from 5 to 90. The number of iterations does not monotonically increase with respect to the network scale although it exhibits an increasing tendency with respect to the network scale. For example, when the level of grid network increases from 20 to 30, the number of iterations decreases from 23 to 21. The relative gaps also are not significantly influenced by the network scale. For example, when the level of grid network increases from 50 to 90, the relative gaps decrease from 4.28 to 3.19.Figure5 shows the computation times relative to the number of nodes of grid network. Figure 6 shows the number of iterations relative to the number of nodes of grid network. As shown in Figure 5, the computational times almost monotonically increase with respect to the number of nodes and links of network, and this indicates that the network scale significantly influences the computational efficiency of the algorithm. This is consistent with the expectations. Increased computational efforts are required for shortest path ork-shortest path searches with increases in the network scale.Figure 5
Computation times relative to the number of nodes.Figure 6
Number of iterations relative to the number of nodes.However, as shown in Figure6, the number of iterations does not monotonically increase with respect to the number of nodes, and this implies that the number of nodes is an important factor that influences the number of iterations, and the number of iterations is also significantly influenced by other initial parameters such as the initial upper bound and lower bound.


## 6.1. A Small Network

A small network is constructed with 5 nodes and 6 links as shown in Figure2. Link travel time tij is a random variable with its mean cij and standard deviation σij. Additionally, each link possesses two resource weights wij1 and wij2. The parameters of each link [cij,σij2,wij1,wij2] are shown in Figure 1.Figure 1
Small network.Figure 2
Updated process of the upper bound, lower bound, and relation gap with respect to the iterations.The initial parameters are set as the OD pair:1-5,ε^=1,λ=0.5, and K=100. The cases of different constraints are computed. Table 1 shows the updated process of Lagrangian multiplier, bound, and gap based on iterations with two upper values of the resource weights: W1=10 and W2=20. Figure 3 shows the updated process of the bound and relative gaps with respect to the iteration. Table 2 shows the updated process of the bound and gap based on iterations with W1=10 and W2=10. Table 3 shows the updated process of the bound and gap based on iterations with W1=20 and W2=20.Table 1
Updated process of the Lagrangian multiplier, bound, and gap with respect to the iterations (W1=10,W2=20).
k μ 1 k μ 2 k μ 3 k U k L k U k - L k ε k 1 100.0 100.0 100.0 60.00 0.00 60.00 100.0  2 61.31 61.48 62.47 52.26 27.20 25.06 47.95  3 63.49 68.23 61.10 48.20 31.12 17.08 35.43  4 53.49 66.15 51.29 47.40 35.49 11.91 25.12  5 51.43 59.30 53.20 43.26 35.98 7.28 16.82  6 30.20 35.36 50.20 42.16 38.32 3.84 9.10  7 34.23 0.00 50.20 40.00 39.13 0.87 2.17  8 32.45 0.00 50.20 40.00 40.00 0.00 0.00Table 2
Updated process of the bound and gap with respect to the iterations (W1=10,W2=10).
k U k L k U k - L k ε k 1 60.00 0.00 60.00 100  2 58.26 29.20 29.06 49.87  3 48.93 36.82 12.11 24.74  4 47.73 45.49 2.24 4.69  5 47.73 47.73 0.00 0.00Table 3
Updated process of the bound and gap with respect to the iterations (W1=20,W2=20).
k U k L k U k - L k ε k 1 60.00 0.00 60.00 100.00  2 52.75 23.20 29.55 56.01  3 45.20 26.12 19.08 42.21  4 43.53 25.40 18.13 41.64  5 38.81 34.33 4.48 11.54  6 38.87 38.87 0.00 0.00Figure 3
Sioux Falls network.As shown in Table1, the exact optimal value is 40, obtained at the 8th iteration. The two resource weights of optimal path 1-3-4-5 are 9 and 12, and this satisfies two constraints. As shown in Figure 2, the upper bounds are gradually close to the lower bounds, and an exact optimal value is finally determined when the relation gap is reduced to 0 in the 8th iteration. As shown in Table 2, an exact optimal value is 47.73, and this is obtained at the 5th iteration. The two resource weights of optimal path 1-2-5 are 7 and 9, and this satisfies two constraints. As shown in Table 3, an exact optimal value is 38.87, and this is obtained at the 6th iteration. The two resource weights of the optimal path 1-2-4-5 are 11 and 12, and this satisfies two constraints. There are differences in the optimal value and optimal path with respect to the constraint conditions.


## 6.2. Sioux Falls Network

Sioux Falls network is a medium-scale transportation network that includes 24 zones, 24 nodes, and 76 links as shown in Figure3 [49].


### 6.2.1. Experimental Design

The experimental design is set as follows. In the Sioux Falls network, link travel timestij((i,j)∈A) are random variables, its mean is cij((i,j)∈A), and its standard deviation is σij((i,j)∈A). In order to simulate its stochastic properties, the expected values and variances are randomly obtained from certain intervals as shown in Table 4. Each link possesses a resource weight wij that is also randomly obtained from a certain interval.Table 4
Intervals of randomly generated link parameters for different links.
c i j σ i j w i j ( i , j ) [ 1,18 ] [ 1,25 ] [ 1,30 ] ( 3,1 ),(3,4),(3,12),(9,5),(9,8),(9,10),(24,13),(24,21),(24,23)  [ 1,25 ] [ 1,25 ] [ 1,20 ] ( 5,4 ),(5,6),(5,9),(10,9),(10,11),(10,15),(10,16),(10,17)  [ 1,20 ] [ 1,28 ] [ 1,25 ] ( 15,10 ),(15,14),(15,19),(15,22),(20,18),(20,19),(20,21),(20,22)  [ 1,18 ] [ 1,20 ] [ 1,19 ] ( 4,3 ),(4,5),(4,11),(8,6),(8,7),(8,9),(8,16),(12,3),(12,11),(12,13),(16,8),(16,10),(16,17),(16,18),  [ 1,15 ] [ 1,30 ] [ 1,20 ] ( 7,8 ),(7,18),(14,11),(14,15),(14,23),(21,20),(21,22),(21,24)  [ 1,20 ] [ 1,20 ] [ 1,20 ] others


### 6.2.2. Computational Results of Different OD Pairs

Table5 shows the computational results of different OD pairs in which the initial parameters are set as ε^=5, W=50, λ=0.2, K=100. The 10 OD pairs are randomly selected from all 24 nodes. As shown in Table 5, all relative gaps are less than the given value ε^=5%, and this implies that the algorithms are terminated within the maximum step K=100. Two relative gaps are zero, and this implies that the accurate solution is finally obtained. All path lengths are less than the constraint value W=50 and satisfy the path length constraint, and this verifies the feasibility and accuracy of the proposed algorithm.Table 5
Computational results of the upper bound, lower bound, and relation gap with respect to the different OD pairs.
OD U L ε Optimal path 1 - 3 41.32 40.45 2.10 1 - 3 - 12 - 11 - 14 - 23 - 24 - 12 - 11 - 14 - 23 - 24  1 - 24 61.94 58.87 4.95 1 - 24  1 - 24 57.31 55.6 2.98 1 - 24  1 - 24 56.43 55.5 1.64 1 - 24  1 - 24 59.24 57.91 2.24 1 - 24  1 - 24 51.66 50.02 3.17 1 - 24  1 - 24 40.60 40.60 0.00 1 - 24  1 - 24 49.13 46.89 4.55 1 - 24  1 - 24 49.08 49.08 0.00 1 - 24  1 - 24 52.14 51.23 1.74 1 - 24


### 6.2.3. Computational Results of Different Constraints

Tables6–8 show the computational results of different constraints with the same OD pair as follows:1-24, 7-19, and 2-23, respectively. The notation NO in Tables 6 and 7 indicates that the algorithm is not successful in obtaining the optimal value when the constrained value W is less than the shortest path weight, a path that satisfies the constraint is absent, and the constrained robust shortest path problem does not include a feasible solution. In Tables 6–8, W=+∞ implies that the constrained robust shortest path problem is reduced to the reliable shortest path problem without constraints. The optimal path varies with respect to the constraint condition that is denoted by constraint value. For example, in Table 6, the optimal paths are different when the constraint values are W=50 and60. Similarly, in Table 7, the optimal paths are also different when the constraint values are W=30 and40. Finally, in Table 8, the optimal paths are also different when the constraint values are W=40,50, and55. Specifically, differences exist in the optimal paths solved by the reliable shortest path problem with and without constraint. For example, in Table 8, the optimal paths are different when the constraint values are W=55 and +∞. Conversely, with increases in the constraint value and path length, the upper bound and lower bound of the optimal path exhibit a decreased trend. For example, in Table 8, the constraint value increases from 40 to 55 and the corresponding path length increases from 36 to 54 although the upper bound decreases from 72.44 to 57.12 and the lower bound decreases from 70.39 to 54.39. The relative gap is not significantly influenced by the constraint condition.Table 6
Computational results of the upper bound, lower bound, and relation gap with respect to the different constraints (1-24).
W U L ε Optimal path 20 NO NO NO NO  50 41.32 40.45 2.10 1 - 3 - 12 - 11 - 14 - 23 - 24  60 37.56 36.12 3.83 + ∞  + ∞ 37.56 37.56 0.00 7 - 19Table 7
Computational results of the upper bound, lower bound, and relation gap with respect to the different constraints (7-19).
W U L ε Optimal path 20 NO NO NO NO  30 51.31 51.48 -0.03 7 - 18 - 16 - 17 - 19  40 33.49 33.49 0.00 7 - 18 - 20 - 19  + ∞ 33.49 33.49 0.00 7 - 18 - 20 - 19Table 8
Computational results of the upper bound, lower bound, and relation gap with respect to the different constraints (2-23).
W U L ε Optimal path 40 72.44 70.39 2.83 2 - 6 - 5 - 4 - 11 - 14 - 23  50 61.94 58.87 4.95 55  55 57.12 54.39 4.77 2 - 1 - 3 - 12 - 11 - 14 - 23  + ∞ 53.20 53.20 0.00 2 - 1 - 3 - 12 - 13 - 24 - 23


### 6.2.4. Computational Results of Different Reliability Coefficients

Table9 shows the computational results of different reliability coefficients in which the initial parameters are set as the OD pair as follows: 14-20, ε^=5, W=50, and K=50. With respect to the first case in Table 9, λ=0 implies that the constrained robust shortest path problem is reduced to the constrained shortest path problem and differences exist between the optimal paths solved by the constrained robust shortest path problem and constrained shortest path problem. The optimal path varies with respect to the reliability coefficient that denotes travelers’ attitudes to path reliability. For example, in Table 9, the optimal paths are different when the reliability coefficient λ=0.2, 0.5, and 0.8.Table 9
Computational results of the upper bound, lower bound, and relation gap with respect to the different reliability coefficients.
λ U L ε Optimal path 0.0 31.31 31.31 0.00 14 - 15 - 19 - 20  0.2 35.12 35.01 0.31 14 - 15 - 19 - 20  0.5 37.25 37.25 0.00 14 - 15 - 22 - 20  0.8 53.45 51.16 4.28 14 - 23 - 22 - 21 - 20


## 6.2.1. Experimental Design

The experimental design is set as follows. In the Sioux Falls network, link travel timestij((i,j)∈A) are random variables, its mean is cij((i,j)∈A), and its standard deviation is σij((i,j)∈A). In order to simulate its stochastic properties, the expected values and variances are randomly obtained from certain intervals as shown in Table 4. Each link possesses a resource weight wij that is also randomly obtained from a certain interval.Table 4
Intervals of randomly generated link parameters for different links.
c i j σ i j w i j ( i , j ) [ 1,18 ] [ 1,25 ] [ 1,30 ] ( 3,1 ),(3,4),(3,12),(9,5),(9,8),(9,10),(24,13),(24,21),(24,23)  [ 1,25 ] [ 1,25 ] [ 1,20 ] ( 5,4 ),(5,6),(5,9),(10,9),(10,11),(10,15),(10,16),(10,17)  [ 1,20 ] [ 1,28 ] [ 1,25 ] ( 15,10 ),(15,14),(15,19),(15,22),(20,18),(20,19),(20,21),(20,22)  [ 1,18 ] [ 1,20 ] [ 1,19 ] ( 4,3 ),(4,5),(4,11),(8,6),(8,7),(8,9),(8,16),(12,3),(12,11),(12,13),(16,8),(16,10),(16,17),(16,18),  [ 1,15 ] [ 1,30 ] [ 1,20 ] ( 7,8 ),(7,18),(14,11),(14,15),(14,23),(21,20),(21,22),(21,24)  [ 1,20 ] [ 1,20 ] [ 1,20 ] others


## 6.2.2. Computational Results of Different OD Pairs

Table5 shows the computational results of different OD pairs in which the initial parameters are set as ε^=5, W=50, λ=0.2, K=100. The 10 OD pairs are randomly selected from all 24 nodes. As shown in Table 5, all relative gaps are less than the given value ε^=5%, and this implies that the algorithms are terminated within the maximum step K=100. Two relative gaps are zero, and this implies that the accurate solution is finally obtained. All path lengths are less than the constraint value W=50 and satisfy the path length constraint, and this verifies the feasibility and accuracy of the proposed algorithm.Table 5
Computational results of the upper bound, lower bound, and relation gap with respect to the different OD pairs.
OD U L ε Optimal path 1 - 3 41.32 40.45 2.10 1 - 3 - 12 - 11 - 14 - 23 - 24 - 12 - 11 - 14 - 23 - 24  1 - 24 61.94 58.87 4.95 1 - 24  1 - 24 57.31 55.6 2.98 1 - 24  1 - 24 56.43 55.5 1.64 1 - 24  1 - 24 59.24 57.91 2.24 1 - 24  1 - 24 51.66 50.02 3.17 1 - 24  1 - 24 40.60 40.60 0.00 1 - 24  1 - 24 49.13 46.89 4.55 1 - 24  1 - 24 49.08 49.08 0.00 1 - 24  1 - 24 52.14 51.23 1.74 1 - 24


## 6.2.3. Computational Results of Different Constraints

Tables6–8 show the computational results of different constraints with the same OD pair as follows:1-24, 7-19, and 2-23, respectively. The notation NO in Tables 6 and 7 indicates that the algorithm is not successful in obtaining the optimal value when the constrained value W is less than the shortest path weight, a path that satisfies the constraint is absent, and the constrained robust shortest path problem does not include a feasible solution. In Tables 6–8, W=+∞ implies that the constrained robust shortest path problem is reduced to the reliable shortest path problem without constraints. The optimal path varies with respect to the constraint condition that is denoted by constraint value. For example, in Table 6, the optimal paths are different when the constraint values are W=50 and60. Similarly, in Table 7, the optimal paths are also different when the constraint values are W=30 and40. Finally, in Table 8, the optimal paths are also different when the constraint values are W=40,50, and55. Specifically, differences exist in the optimal paths solved by the reliable shortest path problem with and without constraint. For example, in Table 8, the optimal paths are different when the constraint values are W=55 and +∞. Conversely, with increases in the constraint value and path length, the upper bound and lower bound of the optimal path exhibit a decreased trend. For example, in Table 8, the constraint value increases from 40 to 55 and the corresponding path length increases from 36 to 54 although the upper bound decreases from 72.44 to 57.12 and the lower bound decreases from 70.39 to 54.39. The relative gap is not significantly influenced by the constraint condition.Table 6
Computational results of the upper bound, lower bound, and relation gap with respect to the different constraints (1-24).
W U L ε Optimal path 20 NO NO NO NO  50 41.32 40.45 2.10 1 - 3 - 12 - 11 - 14 - 23 - 24  60 37.56 36.12 3.83 + ∞  + ∞ 37.56 37.56 0.00 7 - 19Table 7
Computational results of the upper bound, lower bound, and relation gap with respect to the different constraints (7-19).
W U L ε Optimal path 20 NO NO NO NO  30 51.31 51.48 -0.03 7 - 18 - 16 - 17 - 19  40 33.49 33.49 0.00 7 - 18 - 20 - 19  + ∞ 33.49 33.49 0.00 7 - 18 - 20 - 19Table 8
Computational results of the upper bound, lower bound, and relation gap with respect to the different constraints (2-23).
W U L ε Optimal path 40 72.44 70.39 2.83 2 - 6 - 5 - 4 - 11 - 14 - 23  50 61.94 58.87 4.95 55  55 57.12 54.39 4.77 2 - 1 - 3 - 12 - 11 - 14 - 23  + ∞ 53.20 53.20 0.00 2 - 1 - 3 - 12 - 13 - 24 - 23


## 6.2.4. Computational Results of Different Reliability Coefficients

Table9 shows the computational results of different reliability coefficients in which the initial parameters are set as the OD pair as follows: 14-20, ε^=5, W=50, and K=50. With respect to the first case in Table 9, λ=0 implies that the constrained robust shortest path problem is reduced to the constrained shortest path problem and differences exist between the optimal paths solved by the constrained robust shortest path problem and constrained shortest path problem. The optimal path varies with respect to the reliability coefficient that denotes travelers’ attitudes to path reliability. For example, in Table 9, the optimal paths are different when the reliability coefficient λ=0.2, 0.5, and 0.8.Table 9
Computational results of the upper bound, lower bound, and relation gap with respect to the different reliability coefficients.
λ U L ε Optimal path 0.0 31.31 31.31 0.00 14 - 15 - 19 - 20  0.2 35.12 35.01 0.31 14 - 15 - 19 - 20  0.5 37.25 37.25 0.00 14 - 15 - 22 - 20  0.8 53.45 51.16 4.28 14 - 23 - 22 - 21 - 20


## 6.3. Large-Scale Network

### 6.3.1. Experimental Design

The grid network is a scalable network that is typically used to test traffic network equilibrium problem [6]. A grid network of level h contains N=h2 nodes and A=2h(h-1) links as shown in Figure 4: h=4. The experimental design is set as follows. In the grid network, link travel times tij((i,j)∈A) are random variables, its mean is cij((i,j)∈A), and its standard deviation is σij((i,j)∈A). In order to simulate its stochastic properties, all means and standard variances are randomly selected from the interval [0,20]. Each link possesses a resource weight wij that is also randomly selected from the interval [0,20].Figure 4
Grid network.The length constraint of pathW is randomly selected from the interval [0,40×h]. The initial parameters are set as ε^=5, W=50, K=50. The computational times of the algorithm are analyzed under different network scales. All CPU times were measured in seconds with the utility routine CPU time.


### 6.3.2. Computational Time of the Proposed Algorithm

Table10 shows the average computational time of the proposed algorithm for different network scales. The network scale is represented by the level, number of nodes, and links of the grid network. It should be noted that the relative gap ε is the average of 10 computational results with the same OD pairs. The number of iterations k and CPU time is also the average of 10 computational results for the same OD pairs.Table 10
Computational results of different reliability coefficients.
Network scale ε k CPU time (sec.) k N A 5 25 40 0.00 7 4  10 100 180 0.78 11 19  20 400 760 2.41 23 52  30 900 1740 2.23 21 96  50 2500 4900 4.28 41 291  70 4900 9600 4.02 43 485  90 8100 16020 3.19 59 805Evidently, as shown in Table10, the relative gaps are all less than the initial parameters ε^=5. The CPU time significantly increases with respect to the network scale. For example, the CPU time increases from 4 S to 805 S when the level of grid network increases from 5 to 90. The number of iterations does not monotonically increase with respect to the network scale although it exhibits an increasing tendency with respect to the network scale. For example, when the level of grid network increases from 20 to 30, the number of iterations decreases from 23 to 21. The relative gaps also are not significantly influenced by the network scale. For example, when the level of grid network increases from 50 to 90, the relative gaps decrease from 4.28 to 3.19.Figure5 shows the computation times relative to the number of nodes of grid network. Figure 6 shows the number of iterations relative to the number of nodes of grid network. As shown in Figure 5, the computational times almost monotonically increase with respect to the number of nodes and links of network, and this indicates that the network scale significantly influences the computational efficiency of the algorithm. This is consistent with the expectations. Increased computational efforts are required for shortest path ork-shortest path searches with increases in the network scale.Figure 5
Computation times relative to the number of nodes.Figure 6
Number of iterations relative to the number of nodes.However, as shown in Figure6, the number of iterations does not monotonically increase with respect to the number of nodes, and this implies that the number of nodes is an important factor that influences the number of iterations, and the number of iterations is also significantly influenced by other initial parameters such as the initial upper bound and lower bound.


## 6.3.1. Experimental Design

The grid network is a scalable network that is typically used to test traffic network equilibrium problem [6]. A grid network of level h contains N=h2 nodes and A=2h(h-1) links as shown in Figure 4: h=4. The experimental design is set as follows. In the grid network, link travel times tij((i,j)∈A) are random variables, its mean is cij((i,j)∈A), and its standard deviation is σij((i,j)∈A). In order to simulate its stochastic properties, all means and standard variances are randomly selected from the interval [0,20]. Each link possesses a resource weight wij that is also randomly selected from the interval [0,20].Figure 4
Grid network.The length constraint of pathW is randomly selected from the interval [0,40×h]. The initial parameters are set as ε^=5, W=50, K=50. The computational times of the algorithm are analyzed under different network scales. All CPU times were measured in seconds with the utility routine CPU time.


## 6.3.2. Computational Time of the Proposed Algorithm

Table10 shows the average computational time of the proposed algorithm for different network scales. The network scale is represented by the level, number of nodes, and links of the grid network. It should be noted that the relative gap ε is the average of 10 computational results with the same OD pairs. The number of iterations k and CPU time is also the average of 10 computational results for the same OD pairs.Table 10
Computational results of different reliability coefficients.
Network scale ε k CPU time (sec.) k N A 5 25 40 0.00 7 4  10 100 180 0.78 11 19  20 400 760 2.41 23 52  30 900 1740 2.23 21 96  50 2500 4900 4.28 41 291  70 4900 9600 4.02 43 485  90 8100 16020 3.19 59 805Evidently, as shown in Table10, the relative gaps are all less than the initial parameters ε^=5. The CPU time significantly increases with respect to the network scale. For example, the CPU time increases from 4 S to 805 S when the level of grid network increases from 5 to 90. The number of iterations does not monotonically increase with respect to the network scale although it exhibits an increasing tendency with respect to the network scale. For example, when the level of grid network increases from 20 to 30, the number of iterations decreases from 23 to 21. The relative gaps also are not significantly influenced by the network scale. For example, when the level of grid network increases from 50 to 90, the relative gaps decrease from 4.28 to 3.19.Figure5 shows the computation times relative to the number of nodes of grid network. Figure 6 shows the number of iterations relative to the number of nodes of grid network. As shown in Figure 5, the computational times almost monotonically increase with respect to the number of nodes and links of network, and this indicates that the network scale significantly influences the computational efficiency of the algorithm. This is consistent with the expectations. Increased computational efforts are required for shortest path ork-shortest path searches with increases in the network scale.Figure 5
Computation times relative to the number of nodes.Figure 6
Number of iterations relative to the number of nodes.However, as shown in Figure6, the number of iterations does not monotonically increase with respect to the number of nodes, and this implies that the number of nodes is an important factor that influences the number of iterations, and the number of iterations is also significantly influenced by other initial parameters such as the initial upper bound and lower bound.


## 7. Concluding Remarks

The constrained robust shortest path problem is first discussed to obtain the mean-standard deviation shortest path subject to a few resource constraints in a stochastic network in which a reliable objective function and a few resource constraints are simultaneously considered. The constrained shortest path problem is extended from a deterministic network to a stochastic network. It obtainsa priorireliable shortest path subject to resource constraints in the stochastic network. The model is a truer reflection of the actual traffic routing selection than a single model. However, the computational complexity of the model significantly exceeds that of a single model due to the nonlinear and nonadditive nature of the objective function and a few complex constraints.( 1 ) The constrained robust shortest path problem is formulated as a mixed 0-1 convex integer programming model. The objective function is defined as a linear combination of the mean and standard deviation to directly reflect the reliability of travel time. The resource constraint is defined in terms of inequality constraints. The proposed problem is complex due to the nonlinear and nonadditive nature of the objective function, and the inequality constraints that are extremely difficult to calculate in reality.( 2 ) The dual problem of constrained robust shortest path problem is derived based on Lagrangian duality theory. The dual problem is divided into two independent subproblems. The first subproblem corresponds to the classical shortest path problem that is easy to solve. The second subproblem is a nonlinear concave optimization problem, and it is easily solved through minimizing the concave function value at two interval endpoints.( 3 ) An efficient algorithm is constructed to solve a constrained robust shortest path problem. A gradient descent algorithm is developed to solve the dual problem by gradually reducing the gap between the upper and lower bounds. The first subproblem and second subproblem are solved in iterations. The sum of their corresponding objective values is set as the lower bound. Thekth shortest path algorithm is used to obtain the feasible path, and the corresponding objective values are set as the lower bound. When the gap corresponds to the minimum gap or zero, the corresponding upper bound is a near-optimal solution or an accurate solution.( 4 ) An extension of constrained robust shortest path problem that considers multiple constraints is further discussed, and the corresponding dual problem is also derived. An efficient algorithm is introduced to solve this via modifying the proposed algorithm.( 5 ) A numerical test is developed on a small network, medium-scale transportation network, and large-scale network. The numerical results are analyzed. Numerical results indicate the following: (1) the proposed algorithm solves the constrained robust shortest path problem and obtains the near-optimal solution or accurate solution in different scale of networks. This indicates the effectiveness of the proposed algorithm. (2) the computational times of proposed algorithm almost increase linearly with respect to the number of nodes. This indicates the stability of the proposed algorithm. (3) there are differences between the reliable shortest paths obtained under resource unconstraint and constraint, differences also exist between reliable shortest paths obtained under different resource constraints, and the resource constraint significantly influences the choice of reliable shortest path.The constrained robust shortest path problem can be very useful for route guidance systems and logistics optimization problems with complex constraints and reliable-based objectives in stochastic networks. In particular, due to the limitation of battery capacity, the electric vehicle navigation problem is a typical constraint robust shortest path problem. The algorithm can be applied to real road network as the core of in-vehicle route guidance systems in the future.On the basis of the algorithm proposed in this paper, some extensions can be envisaged. First, in this study link travel times are assumed to follow a spatial independent distribution. Previous studies have found that travel times could be better approximated by correlated link travel times [6]. In future study correlations between link travel times can be incorporated either through statistical model or through traffic network simulation to capture a traffic network in a more realistic way. Second, in this study link travel times are assumed to be time independent; dynamic characteristics among traffic network are not incorporated in the model. Actually link travel times are not only random but also dynamic in realistic traffic network [50]. The constrained robust shortest path problem should be extended from static network to dynamic network. Third, the constrained robust shortest path problem was formulated based on the mean-standard deviation shortest path concept. However, other models of reliable shortest path have been developed in different application contexts, such as α-reliable path [17]. The extension of the constrained robust shortest path problem to other reliable shortest path models will be another topic for further studies.



---

*Source: 3987278-2019-06-27.xml*