# Estimating the Capacity of Urban Transportation Networks with an Improved Sensitivity Based Method

Authors: Muqing Du; Xiaowei Jiang; Lin Cheng
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

The throughput of a given transportation network is always of interest to the traffic administrative department, so as to evaluate the benefit of the transportation construction or expansion project before its implementation. The model of the transportation network capacity formulated as a mathematic programming with equilibrium constraint (MPEC) well defines this problem. For practical applications, a modified sensitivity analysis based (SAB) method is developed to estimate the solution of this bilevel model. The high-efficient origin-based (OB) algorithm is extended for the precise solution of the combined model which is integrated in the network capacity model. The sensitivity analysis approach is also modified to simplify the inversion of the Jacobian matrix in large-scale problems. The solution produced in every iteration of SAB is restrained to be feasible to guarantee the success of the heuristic search. From the numerical experiments, the accuracy of the derivatives for the linear approximation could significantly affect the converging of the SAB method. The results also show that the proposed method could obtain good suboptimal solutions from different starting points in the test examples.


---

## Body

## 1. Introduction

As the rapid growth of urbanization, the population and economy of most cities in the developing countries or regions are going through significant changes. The fast development of transportation infrastructures in these areas gives rise to a quick change on the travel behaviors. Along with the city’s expansion, an inevitable thing is to design new transportation system which is capable of meeting the future development in land-using and population growth. However, under fast-changing of travel demand, the conventional forecasting methods (e.g., “four-step”) could not provide a straightforward evaluation to the new designed transportation network such as on whether the new network is suitable for the travel demand in the target year or how many trips can be accommodated. Thus, to help the government to make decision on the expansion or construction project or preassessment of social benefits, the throughput, or “capacity,” of the given transportation network is of practical meaning to be estimated before planning implementation.Capacity is a commonly used property to represent the maximum flows that can pass through the link or node in transportation system. In an attempt to address the question of what is the maximum attainable throughput of the given network, the concept of capacity is employed as an important measurement for transportation system evaluation. It is able to reflect how much traffic demand can be accommodated by a given transportation system. Thus, efficient policy for land use or traffic restraint and growth can be established in advance. According to the conventional network flow theory, the capacity problem is stated to find the maximum flows that can be sent from a specified source node to another specified sink node without exceeding the capacity of any link. This well-known problem is extended to the multicommodity and is widely used in freight transport. However, when modeling the capacity of urban transportation network, the problem becomes quite complex. Noted by Yang et al. [1], travelers in urban transportation network can choose their routes and their trip costs increase with increasing flow as a result of congestion. Besides, multiple origin and destination (O-D) pairs exist and the flows between distinct O-D pairs cannot be exchanged in passenger transportation system. These differences make the modeling of transportation network capacity complex, and the intriguing problem is also hard to solve.The most popular formulation of the transportation network capacity is the bilevel programming model, which maximizes the traffic flows under the equilibrium constraint. Wong and Yang [2] first incorporated thereserve capacity concept into a traffic signal control network. This concept was widely extended in the study of signal controlled networks [3–5]. The reserve capacity is defined as the largest multiplier applied to a given O-D demand matrix that can be allocated to a network, so the solution is significantly affected by the predetermined O-D matrix. However, it is unrealistic to assume that all O-D flows increase in the same rate, especially for the areas under rapid changing. If the predetermined distribution proportion is far from the future tendency, the solution will be of little use. Consequently, in order to reflect the differences of the future development of each urban subarea, Yang et al. [1] considered that the new increased O-D demand pattern should be variable both in level and in distribution, while the distribution of the current trips would be relatively fixed. Thereby they introduced the equilibrium trip distribution/assignment model with variable destination costs (ETDA-VDC) [6] to capture this characteristic for network capacity estimation. Based on this model, Kasikitwiwat and Chen [7] proposed the concepts and models of theultimate andpractical capacity. The former is used for the new network without any current flow, while the latter is the same as Yang’s model. Then, Chen and Kasikitwiwat [8] used the practical network capacity model to describe the limited flexibility of transportation networks. According to the literatures, the concept of the practical capacity model is more fully functional and preferred, as it takes both the current demand pattern and the variability of future growth into consideration.In order to solve the various bilevel capacity models, the SAB algorithm is generally employed. This SAB algorithm for bilevel programing was first presented in Friesz et al. [9]. It is heuristic and depends absolutely on the derivative information produced by the sensitivity analysis of the lower-level problem [10]. Benefiting from the rich achievements in the study of the sensitivity analysis for equilibrium models [11–18], SAB algorithm has been widely utilized in the optimization problems of equilibrium network flows, such as traffic signal control [3, 5, 19, 20] and network design [21–23], as well as network capacity [1, 3–5, 24]. But due to the difficulty of the sensitivity analysis for ETDA-VDC model, Yang et al. [1] used an iterative estimation-assignment (IEA) algorithm [11] instead to solve the transportation network capacity problem. Later, Kasikitwiwat and Chen [7] and Chen and Kasikitwiwat [8] selected using a genetic algorithm to solve the problem in a very small network. However, since the complexity of the network capacity problem, the global optimization algorithms (e.g., genetic algorithm or simulated annealing) can hardly find the exact solution to the capacity problem in larger networks, and the computation time could be intolerably long. By contrast, SAB algorithm has the property of fast convergence which makes the computation terminate at a local optimum within a considerable time. Nevertheless, the calculation issue of the matrix inversion still limits the applications of the SAB method. To address this problem, we developed an effective method by simplifying the matrix inversion in the sensitivity analysis approach, which will take much less memory space, so the capacity of the real transportation networks could be estimated.In this study, Yang’s formulation [1] for the transportation network capacity model is employed to describe the practical capacity of the urban road system. In an attempt to estimate the capacity of the real road networks, a series of improvements are taken to the SAB method to make the heuristic search successfully converge to a relatively better suboptimal solution. Firstly, the OB algorithm [25] is modified for the solution of the lower-level ETDA-VDC model. Then, the restriction sensitivity analysis approach for the ETDA-VDC model [18] is employed in and improved on the expressions so as to deal with the large-scale problems. Besides, the solution update strategy is modified on the step-size adaption, which ensures the entire heuristic search to converge to a local optimum. Finally, numerical experiments are implemented to show the efficiency and capability of the proposed SAB method.


## 2. Road Network Capacity Model

It should be noted that the boldface type of the Notation section represents the corresponding column vectors in the remainder of this paper.


### 2.1. Model Formulation

Conversional methods, like the reserve capacity model, evaluate the capacity of transportation networks by assuming that the travel demand increases with a determined distribution proportion, which is usually far from the regularity and underestimates the results. In order to evaluate zonal development potential and equilibrium network capacity more appropriately, Oppenheim’s definitions on the behaviors of the existing travel demand and the additional demand are introduced as follows.(i)
The existing demand, denoted by epq, has predetermined origins and destinations. The pattern of the existing demand is formed during the past long term, so its distribution is going to be relatively stable and can be regarded as fixed. The existing demand only changes routes to optimize the travel cost.(ii)
The additional demand, dpq, is variable. The new generated demands from residential area can decide their daily travels without the constraints of either destination or route choices. But the behavior of the additional demand still follows the rule that selects the destinations which maximize the “utility” of the trips.The utility could include the destination attractiveness, the cost along traveling route, and other factors. The attractiveness of destination is determined by the congestion at destination and the expenses for the activity in that area. In network capacity model, the utility from originp to destination q is formulated as Upq=-(Mq+τpq), in which τpq is the travel cost from p to q, and Mq denotes the destination cost which could be a decrease function of the total additional trip attraction, Dq, at destination q. Besides, the destination choices of the travelers at each origin p are assumed to have certain randomness. Thus, the conditional probability that an individual will choose destination p is derived by using the standard logit function, so the O-D travel demand is conducted by(1)dpq=Opexp⁡-θτpq+Mq∑k∈Zexp⁡-θτpk+Mq,∀p∈Z,q∈Zp.Thus, with the objective to maximize the additional demand under the above travel behavior regularity and certain physical constraints, the typical road network capacity model is formulated as the following bilevel programming problem.Upper-level problemis as follows:(2)Max⁡O∑p∈ZOp,(3)s.t.vaO≤Ca,∀a∈A,(4)Op=∑q∈ZpdpqO≤Opmax⁡-O-p,∀p∈Z,(5)Dq=∑p∈ZqdpqO≤Dqmax⁡-D-q,∀q∈Z,(6)Op≥0,∀p∈Z,where dpq(O) and va(O) are obtained by solving the ETDA-VDC problem in lower-level problem.Lower-level problem, ETDA-VDC model,is as follows:(7)Min⁡(f,h,d)∑a∫0vataxdx+1θ∑p∈Z∑q∈Zpdpq(ln⁡dpq-1)vvvvvvvvvv+∑q∈Z∫0DqMq(y)dy,(8)s.t.∑q∈Zpdpq=Op,∀p∈Z,(9)∑r∈Rpqhrpq=epq,∀p∈Z,q∈Zp,(10)∑r∈Rpqfrpq=dpq,∀p∈Z,q∈Zp,(11)va=∑p∑q∑rfrpq+hrpqδa,rpq,∀a∈Am,(12)dpq≥0,∀p∈Z,q∈Zpm,(13)frpq≥0,∀p∈Z,q∈Zp,r∈Rpq,(14)hrpq≥0,∀p∈Z,q∈Zp,r∈Rpq.The upper-level problem defines a maximal trip production model. The objective is to maximize the summation of the additional trip production at origins. Equation (3) represents that the traffic flow on every link should not exceed its capacity. Constraints (4) and (5) are the limitation of the zonal trip production and attraction. They mean the number of trips generated and attracted at each traffic zone should be limited by some upper bounds, namely, Opmax⁡ and Dqmax⁡, respectively.The lower-level problem is the ETDA-VDC model. The objective function (7) indicates the choice behavior of both the existing and additional travel demand. Constraint (9) shows that the amount of the existing flows is fixed for each O-D, while constraints (8) and (10) show that the additional flows are only restrained at the origin productions. The relationship between the link flow and route flow is represented in (11). All the variables must be nonnegative, that is, constraints (12)–(14). The lower-level problem is a combined distribution and assignment model.This bilevel model was first presented in work by Yang et al. [1]. Because of the advantages on the formulation of the travel demand growth, it was continually used in later researches as a typical model for the road network capacity concept. The remaining part of this study focuses on the solution of this model in the real-sized road networks.


## 2.1. Model Formulation

Conversional methods, like the reserve capacity model, evaluate the capacity of transportation networks by assuming that the travel demand increases with a determined distribution proportion, which is usually far from the regularity and underestimates the results. In order to evaluate zonal development potential and equilibrium network capacity more appropriately, Oppenheim’s definitions on the behaviors of the existing travel demand and the additional demand are introduced as follows.(i)
The existing demand, denoted by epq, has predetermined origins and destinations. The pattern of the existing demand is formed during the past long term, so its distribution is going to be relatively stable and can be regarded as fixed. The existing demand only changes routes to optimize the travel cost.(ii)
The additional demand, dpq, is variable. The new generated demands from residential area can decide their daily travels without the constraints of either destination or route choices. But the behavior of the additional demand still follows the rule that selects the destinations which maximize the “utility” of the trips.The utility could include the destination attractiveness, the cost along traveling route, and other factors. The attractiveness of destination is determined by the congestion at destination and the expenses for the activity in that area. In network capacity model, the utility from originp to destination q is formulated as Upq=-(Mq+τpq), in which τpq is the travel cost from p to q, and Mq denotes the destination cost which could be a decrease function of the total additional trip attraction, Dq, at destination q. Besides, the destination choices of the travelers at each origin p are assumed to have certain randomness. Thus, the conditional probability that an individual will choose destination p is derived by using the standard logit function, so the O-D travel demand is conducted by(1)dpq=Opexp⁡-θτpq+Mq∑k∈Zexp⁡-θτpk+Mq,∀p∈Z,q∈Zp.Thus, with the objective to maximize the additional demand under the above travel behavior regularity and certain physical constraints, the typical road network capacity model is formulated as the following bilevel programming problem.Upper-level problemis as follows:(2)Max⁡O∑p∈ZOp,(3)s.t.vaO≤Ca,∀a∈A,(4)Op=∑q∈ZpdpqO≤Opmax⁡-O-p,∀p∈Z,(5)Dq=∑p∈ZqdpqO≤Dqmax⁡-D-q,∀q∈Z,(6)Op≥0,∀p∈Z,where dpq(O) and va(O) are obtained by solving the ETDA-VDC problem in lower-level problem.Lower-level problem, ETDA-VDC model,is as follows:(7)Min⁡(f,h,d)∑a∫0vataxdx+1θ∑p∈Z∑q∈Zpdpq(ln⁡dpq-1)vvvvvvvvvv+∑q∈Z∫0DqMq(y)dy,(8)s.t.∑q∈Zpdpq=Op,∀p∈Z,(9)∑r∈Rpqhrpq=epq,∀p∈Z,q∈Zp,(10)∑r∈Rpqfrpq=dpq,∀p∈Z,q∈Zp,(11)va=∑p∑q∑rfrpq+hrpqδa,rpq,∀a∈Am,(12)dpq≥0,∀p∈Z,q∈Zpm,(13)frpq≥0,∀p∈Z,q∈Zp,r∈Rpq,(14)hrpq≥0,∀p∈Z,q∈Zp,r∈Rpq.The upper-level problem defines a maximal trip production model. The objective is to maximize the summation of the additional trip production at origins. Equation (3) represents that the traffic flow on every link should not exceed its capacity. Constraints (4) and (5) are the limitation of the zonal trip production and attraction. They mean the number of trips generated and attracted at each traffic zone should be limited by some upper bounds, namely, Opmax⁡ and Dqmax⁡, respectively.The lower-level problem is the ETDA-VDC model. The objective function (7) indicates the choice behavior of both the existing and additional travel demand. Constraint (9) shows that the amount of the existing flows is fixed for each O-D, while constraints (8) and (10) show that the additional flows are only restrained at the origin productions. The relationship between the link flow and route flow is represented in (11). All the variables must be nonnegative, that is, constraints (12)–(14). The lower-level problem is a combined distribution and assignment model.This bilevel model was first presented in work by Yang et al. [1]. Because of the advantages on the formulation of the travel demand growth, it was continually used in later researches as a typical model for the road network capacity concept. The remaining part of this study focuses on the solution of this model in the real-sized road networks.


## 3. Sensitivity Based Heuristic Algorithm

This section presents an improved version of the SAB algorithm for the solution of the road network capacity model. To overcome the drawbacks of the conventional SAB algorithm [9] that cannot be applied to any real-sized network for capacity estimations, the following improvements are carried out.(i)
The lower-level ETDA-VDC model is fast solved by a modified OB algorithm to produce a high accurate solution.(ii)
The rectified sensitivity analysis method for the ETDA-VDC model is simplified on the calculation of the matrix inverse to be applicable for large-scale problems.(iii)
The solution update of the heuristic search is improved by step-size adaption in order to ensure that the SAB algorithm can converge to a local optimum.Correspondingly, a series of techniques is proposed in this section, so our modified SAB algorithm will be capable to solve the bilevel road capacity model efficiently.


### 3.1. Origin-Based Algorithm for the ETDA-VDC Model

In the standard SAB search, the lower-level ETDA-VDC model should be solved in every iteration to conduct an equilibrium traffic flow pattern, namely, the solution to the lower-level problem. According to recent researches on the traffic assignment problems, the OB algorithm is demonstrated to be one of the state-of-the-art algorithms [26]. In addition, from the results of the OB algorithm, the set of all equilibrated routes can be easily extracted, which will be a precondition for the restriction sensitivity analysis approach in next step of SAB algorithm.The OB algorithm uses theorigin-based approach proportions, αp={αap}, as its main solution variables, where αap represents the proportion of flow that comes from origin p through link a. The approach proportions are updated by shifting flows within the restricting subnetwork, Ap. The details of the OB algorithm can be referred to in Bar-Gera [26]. Using the approach proportions, the OB algorithm is able to store the route flow information with significantly less memory than the route-based algorithms and achieve high-accuracy solutions compared to the link-based algorithms.In the process of the OB algorithm, it starts with trees of minimum cost routes from origins as restricting networks. Then steps of updating restricting network and approach proportions are implemented for each origin separately. An “inner loop” is performed subsequently to accelerate convergence, in which the origin-based link flows are updated while keeping the restricting subnetworks fixed. The inner loop is useful, because the restricting subnetworks tend to stabilize fairly quickly but updating the restricting subnetworks requires more computational effort [26, 27].The OB algorithm for the combined trip distribution and assignment problems requires an additional step to update the O-D flows while keeping the route proportions fixed. Taking Evans’ algorithm as a reference, we modified the original OB algorithm for the solution of the ETDA-VDC model, in which the step size of the O-D flow update is obtained by solving a one-dimension search problem. The proposed algorithm is summarized as shown in Algorithm1.Algorithm 1:
Initialization:
Setn∶=0, determine the initial link cost t0 and destination cost M0 by setting x0=0 and d0=0, respectively.
For each originp∈Z do
(i) Find the tree of minimum cost routes rooting fromp. Let Ap be the minimum route tree. Denote
the minimum route cost to destinationq∈Zp by τpq0, and choose a minimum cost route from pto q.
(ii) Compute the initial variable O-D demands by
dpq0=Opexp⁡[-θ(τpq0+Mq0)]∑kexp⁡[-θ(τpk0+Mk0)],∀q∈Zp
(iii) Forq∈Zp, assign the entire O-D demand epq and dpq0 to the minimum cost route r from p to q, and obtain initial
link flowxap0.
(iv) Update the link costs using the initial link flows.
(v) Initialize the origin-based approach proportionsαap0.
Main Loop:
Given the current variable O-D demanddn obtained in (n-1)th-iteration:
Forn=1 to number of main iterations (IMain)
for eachp in Z do
Update restricting networkAp
Update origin-based approach proportionsαp
end for
Inner Loop:
form=1 to number of inner iterations (Iinner)
for eachp in Z do
Update origin-based approach proportionsαp
end for
end for
UpdateO-Dflows, retain origin-based approach proportions:
Given the origin-based approach proportionsαpn, link flows vn and link costs tn obtained in the steps
above:
(i) For eachq∈Z, compute the destination cost Mqn by O-D demand dpqn and epq.
(ii) Find the set of auxiliary trip demandsd^pqn by solving the following logit distribution model:
d^pqn=Opexp⁡[-θ(τpqn+Mqn)]∑kexp⁡[-θ(τpkn+Mkn)]∀p∈Z,q∈Zp
(iii) Calculate the auxiliary traffic flowv^an on each link a with the approach proportions αapn
(iv) Let(vλ(n),dλ(n))=(1-λ)(vn,dn)+λ(v^n,d^n) and solve the one-dimensional search problem defined
as follows to obtain the step sizeλ*∈[0,1]
min0≤λ≤1⁡F(λ)=∑a∈A∫0vaλ(n)ta(x)dx+(1/θ)∑p∈Z∑q∈Zpdpqλ(n)(ln⁡dpqλ(n)-1)+∑q∈Z∫0Dqλ(n)Mq(y)dy
(v) Set(vn+1,dn+1)∶=(vλ(n),dλ(n)) and check for convergence. Terminate if the convergence
criterion is satisfied; otherwise, update total link flows and link costs, setn∶=n+1 and start
a new iteration of the main loop.
End forThe performance of Algorithm1 was testified on different test networks. The characteristics of these networks are obtained from http://www.bgu.ac.il/~bargera/tntp/. The converging processes on two well-known networks are demonstrated in Figure 1, where the relative gap is defined as RG=|Fn-Fn-1|/Fn (Fn is the value of objective function (7) at nth main loop). It shows that the converging of the improved OB algorithm (IOBA) is faster and more stable than Evans’ algorithm (EA) or the original OB algorithm (OBA) steps in Bar-Gera and Boyce [25].Figure 1
Relative gap versus CPU time.
(a)
(b)


### 3.2. Sensitivity Analysis of the ETDA-VDC Model

On the basis of an equilibrium solution of the lower-level model, the derivatives of the lower-level decision variables with respect to the upper-level ones should be produced for linearly approximating the whole bilevel model in next step. The derivative expressions can be obtained by employing the restriction sensitivity analysis approach for the ETDA-VDC model [18]. The restriction approach reduces the original network to a restricting one in which the nonuniqueness difficulty will be overcome, and thus the full derivatives of the link flows in the primal problem can be solved. The necessary results from the sensitivity analysis of the ETDA-VDC model are presented in this section without proof.In the restricting problem the derivatives of the model solutions with respect to the input parameters are derived from(15)∇Ox(ε)=Jx-1-JO,where Jx is an invertible Jacobian matrix with respect to the solution variables x=(fBT,hBT,dT,λT,μT,uT)T, in which fB and hB are the basic variables of the additional and existing route flows, respectively. JO is the Jacobian matrix with respect to the additional zonal productions, O. Here, λ,μ,u are the Lagrange multipliers associated with the conservation constraints (8)–(10). O is referred to as an input parameter for the ETDA-VDC model. It is in terms of the perturbation ε. The expressions of Jx and JO are as follows:(16)Jx=∇f2L∇f,hLOOO-ΛfT∇h,fL∇h2LOO-ΛhTOOO∇d2L-ΦTOIOO-ΦOOOO-ΛhOOOO-ΛfOIOOO,(17)JO=OOOIOOT,where ∇f2L=ΔfT∇vt(v,ε)Δh, ∇h,fL=ΔhT∇vt(v,ε)Δf, ∇h2L=ΔhT∇vt(v,ε)Δh, ∇f,hL=ΔfT∇vt(v,ε)Δh, and ∇d2L=(1/θ)diag(1/dpq)+ΨT∇dM; diag(1/dpq) is a diagonal matrix with 1/dpq as its diagonal elements. The superscript “T” represents the transposed matrix. “I” indicates the identity matrix. The link cost function, t(·), should be strongly monotonically defined in v, which guarantees the uniqueness of the equilibrium solution and is also a precondition for the sensitivity analysis approach. Utilizing the restriction approach, the derivative, ∇Ox(ε), is first conducted in the restricting problem and then ∇Ov, which is defined in the original one, can be obtained by(18)∇Ov=Δf∇OfB+Δh∇OhB.Equation (18) indicates that the variations on the link flows can be just represented by the changes on the basic route flows, which is the rationale of the restriction sensitivity analysis approach.From (15), the inverse of the matrix Jx is required to be worked out. However, the calculation of the inverse matrix is costly for both computation time and storage space in numerical computation for large-scale problems. The dimension of Jx (from (16)) could be about five times of the number of O-D pairs. So the elements of Jx could totally take approximately O(25n4) space, where n is the amount of the traffic zones, and thus even a medium network with hundreds of origins will conduct a Jacobian matrix with over one hundred thousand dimensions. For example, in a small network with about 30 traffic zones (900 O-D pairs), nearly 400 MB RAM should be expended on the storage of the Jacobian matrix (if each element uses a double type). Consequently, the conventional sensitivity analysis approach for the combined model is very difficult to be utilized in the large-scale problems. Besides, the direct calculation for the inverse of a large matrix is pretty inefficient and inaccurate in practice.To avoid the defects of deriving the inverse matrix directly, the usual way to solve the following linear equation as an alternative. Equation (15) which produces the sensitivity results can be rewritten as(19)Jx·∇Ox=-JO.Then, ∇Ox is derived by solving the following series of linear equations:(20)bi=Jx·zi,where bi is the ith column vector in matrix of [-JO] and zi is the ith column vector in ∇Ox. Let ei be a unit vector with one in the ith position and zeros elsewhere. The length of ei is equal to the column number of JO. Thus, (21)bi=000-ei00,-JO=…,bi,…,where 0 in bold is zero vectors associate with the zero block in JO.To solve (20), Jx is premultiplied by the matrices K1 and K2 in sequence and postmultiplied by K1T and K2T simultaneously, which is equivalent to making the elementary transformation of matrix Jx. Consider(22)K1=IIIIΦII,K2=IIIII-∇d2L-1I.Thus,(23)K2K1JxK1TK2T=∇f2L∇f,hLO-ΛfTΦTO-ΛfT∇h,fL∇h2LOO-ΛhTOOO∇d2LOOO-ΦΛfOOOOOO-ΛhOOOO-ΛfOOOO-∇d2L-1.By using the above operations on (20), we obtain(24)K2K1JxK1TK2TK2T-1K1T-1zi=K2K1bi.Obviously, (25)y=K2T-1K1T-1zi,K2K1bi=000-ei00.Therefore, it is equivalent to solving the equation in terms ofy and then deriving zi by(26)zi=K1TK2Ty=y1y2y3-∇d2L-1·y6y4y5y6+ΦT·y4.Since it is only concerned with the first three subvectors ofzi in the sensitivity result, we just need to find the value of y1, y2, y3, and y6. From (24), the following group of equations can be obtained:(27)∇f2L·y1+∇f,hL·y2-ΛfTΦT·y4-ΛfT·y6=0,(28)∇h,fL·y1+∇h2L·y2-ΛhT·y5=0,(29)∇d2L·y3=0,(30)-ΦΛf·y1=-ei,(31)-Λh·y2=0,(32)-Λf·y1-∇d2L-1·y6=0.The above equation systems can be further simplified depending on whether the incidence matrix Λh is square.(1)
IfΛh is square, consider the following.Thus,Λh will be constructed as an identity matrix; the value of y2 and y3 can be computed from (29) and (31), and thereby y2=y3=0. From (32),(33)y6=-∇d2L·Λf·y1.Substituting (33) into (27), we obtain(34)∇f2L+ΛfT∇d2LΛfy1-ΛfTΦT·y4=0.Combined with (30) we get(35)∇f2L+ΛfT∇d2LΛf-ΛfTΦT-ΦΛfOy1y4=0-ei,where(36)∇f2L+ΛfT∇d2LΛf=ΔfTΛfT∇vt(·)∇d2LΔfΛf.Because the matrix ∇vt(·)∇d2L is positive definite, the columns of ΔfΛf are leaner independent from the restriction approach. It can be easily conducted that matrix (36) is invertible. Let(37)∇f2L+ΛfT∇d2LΛf-ΛfTΦT-ΦΛfO=ABTBO.Solving (35), we can get(38)y1=A-1BTBA-1BT-1·ei.To this extent, the only effort to produce the desired derivative results is to calculatey1, in which the inverse of matrices A and BA-1BT should be derived. The dimension of matrix A is equal to the numbers of the basic routes used by the additional travel demand, which is a little more than the number of O-D pairs from our observation in computational experiments. The dimension of matrix BA-1BT is the same as the number of origins and is not large.Consequently, by utilizing (33) and y2=y3=0, the first three subvectors of zi can be obtained, which are separately corresponding to the derivative results, ∇OfB, ∇OhB, and ∇Od, in ∇Ox. The expressions are(39)∇OfB=A-1BTBA-1BT-1,∇OhB=0,∇Od=Λf·∇OfB.Thus, (40)∇Ov=Δf∇OfB.(2)
IfΛh is not square, consider the following.In this case,y2=∇OhB will not be equal to zero. For further simplification, we rewrite Λh as(41)Λh=Λh′I.Let R′ denote the set of routes associated with Λh′ and the flows on these routes are represented by h′. Repeating the derivation from (33) to (35), we can get(42)∇f2L+ΛfT∇d2LΛf∇f,hL′-ΛfTΦTO∇h,fL′∇h2L′OΛh′T-ΦΛfOOOOΛh′TOO×y1y2′y4y5′=00-ei0,where the superscript “′” indicates the matrices and variables corresponding to the routes in R′. Considering that the set of R′ is fairly small for network capacity problems, the equation system is slightly larger than (35), which is also easy to be solved. The derivative of the existing route flows is computed by y2′. Consider(43)∇OfB=y1,(44)∇OhB=y2′,(45)∇Od=Λf·∇OfB.And thus, (46)∇Ov=Δfy1+Δh′y2′.From the above result, it can be noted that the primal problem to calculate the inverse ofthe entire matrixJx is replaced by solving a small system of linear equations.


### 3.3. Restriction of the Equilibrium Network

When the restriction sensitivity analysis approach was proposed [10], the restricting problem was conducted by finding an extreme point of the feasible region of the equilibrium route flow. An equivalent linear programming was provided for solution. However, this method has been verified usually failed, since the linear programming may find a degenerated solution which causes the fact that the reduced restricting problem cannot reserve the essential information and conducts a wrong result [16]. A rectified method is to find a maximum set of the linear independent columns in the coefficient matrix of the following equation system [18]:(47)Δf0Δh0Λf0OOΛh0f0h0=v*d*e.The above equation system defines the feasible region of the equilbirum route flow under unique euqilibrium solution of v* and d*. Note that the notations are distinguishing from the ones in (16). f0and h0 correspond to the routes that carry positive flows, which are the equilibrated (or minimum-cost) routes according to Assumption 2 in Du et al. [18], and can be obtained from the results of the OB algorithm. The superscript “0” is associated with the variables in a reduced problem, which only consists of the equilibrated routes.From the coefficient matrix of (47), the maximum set of the linear independent columns corresponds to the set of the equilibrated and linear independent routes (ELI, denoted by RB) and can be found using the row echelon form of the coefficient matrix. We can utilize the blocks Λf and Λh which are incidence matrices with one in each column (because one route can only serve one O-D). Since there may be more than one equilibrated route between an O-D pair, Λf and Λh are generally not column full rank. Thus, the coefficient matrix can be rewritten as (48)Λ-fOΛ~fOOΛ-hOΛ~hΔ-fΔ-hΔ~fΔ~h=IΛ~fOIOΛ~hΔ-fΔ-hΔ~fΔ~h⟶IOOOΔ~f-Δ-fΛ~fΔ~h-Δ-hΛ~h,where Λ-f and Λ-h are square matrices, which consist of exactly one equilibrated route for each O-D, so they are equal to identity matrices. Thus, the question is reduced to find and eliminate the linear dependent columns in a submatrix. Consider(49)Δ~f-Δ-fΛ~fΔ~h-Δ-hΛ~h.The columns of the above matrix are only related to the alternative routes between every O-D. Therefore, the set RB will include the routes corresponding to Λ-f, Λ-h, and the maximum set of the independent columns in (49).


### 3.4. Solution of the Maximum Trip Production Problem

When the sensitivity results are derived, the SAB algorithm will use this information to represent the implicit relationships betweenv, d, and O. Thus, the bilevel problem is first-order approximated at the given point, O*=f-1(v*,d*). Let v(O*) and d(O*), respectively, denote the solutions to the lower-level model at O*. The relationship can be represented by using the Taylor expansion:(50)vO≈vO*+∇Ov·O-O*,dO≈dO*+∇Od·O-O*,where the derivatives ∇Ov and ∇Od are obtained from the sensitivity analysis of the ETDA-VDC model. Therefore, the bilevel problem can be reformulated as(51)Max⁡∑p∈ZOps.t.∇Ov·O≤C-v*+∇Ov·O*Φ·∇Od·O≤Omax⁡-O--O*+Φ·∇Od·O*Ψ·∇Od·O≤Dmax⁡-D--D*+Ψ·∇Od·O*.The solution of this linear programming can be easily derived using the simplex method. However, because the linear programming problem is just locally approximated, the new solution might be infeasible to the original problem. At an infeasible point of O* the lower-level trip distribution and assignment results may not satisfy the upper-level constraints. In extreme cases, for example, some link flow va* may be much greater than its capacity Ca, which could cause the capacity constraints to fail for any nonnegative O. In consequence, the new linear approximation will have no solution.In consideration of this flaw, we modified the conventional SAB algorithm by restraining the solution to be always located within the feasible region. A trip distribution and assignment step is implemented at the solution of the above linear programming,O*, and then the results are checked with the capacity conditions. If any capacity constraint is violated, the new solution should be updated by a convex combination of the solution from last iteration, O(n), and O*(n). Therefore, let O^ be the solution to the approximate linear problem above; the maximum step size λ=2-k(k=0,1,2,…) is chosen which ensures the capacity constraint is satisfied at any link or traffic zone. The new solution is updated by(52)O(n+1)∶=O(n)+2-λO^n-On.


### 3.5. A Heuristic Solution Process

In summary, the modified SAB algorithm involves an iterative process between the upper-level and lower-level problems and can be summarized as follows.Step
0 (initialization). Determine an initial value of trip production pattern O(0). Set n∶=0.Step  1 (solving lower-level problem). Solve the ETDA-VDC model in lower-level for the given O(n), by which the equilibrium link flows v(n) and O-D demand d(n) are obtained. The set of all equilibrated routes can be derived at the equilibrium point.Step
2 (sensitivity analysis). Calculate the partial derivatives, ∇Ov and ∇Od, using the sensitivity method for the ETDA-VDC model.Step
3 (local linear approximation). Formulate local linear approximations of the upper-level capacity constraints using the derivative information and solve the approximate linear programming problem to produce an auxiliary trip production O^(n).Step
4 (updating solution). Let O(n+1)∶=O(n)+2-λ(O^(n)-O(n)), where λ=0,1,2,…, until the capacity constraints are satisfied by solving the ETDA-VDC problem with O(n+1). Set n:=n+1.Step  5 (convergence criterion). If |Op(n+1)-Op(n)|≤κ for all p∈Z then stop, where κ is a predetermined tolerance. Otherwise, return to Step 1.The above process is a modification of the conventional SAB approach. The flowchart is shown in Figure2. It should be noted that without the step-size adaption in Step 4 the SAB method only works in very small networks but could fail in larger examples. The set of all equilibrated routes is used to apply the restriction approach for sensitivity analysis. In this study, the modified OB algorithm for the ETDA-VDC problem is used for this purpose. In addition, the derivative results of the sensitivity analysis can give a precise local linear approximation to the upper-level capacity constraints, which is very important to make the heuristic search converge to a considerable good solution.Figure 2
Flowchart of sensitivity analysis based algorithm for network capacity model.In addition, given the nonconvexity of the bilevel problem, the SAB algorithm will converge to a local optimal point [2]. However, as shown in Yang et al. [28], for MPEC models, if the upper-level objective function is a linear function of its decision variables, the heuristic algorithm can at least find a noninferior optimal solution. So the SAB algorithm is supposed to be able to find a satisfying solution for the road network capacity problem in practice.


## 3.1. Origin-Based Algorithm for the ETDA-VDC Model

In the standard SAB search, the lower-level ETDA-VDC model should be solved in every iteration to conduct an equilibrium traffic flow pattern, namely, the solution to the lower-level problem. According to recent researches on the traffic assignment problems, the OB algorithm is demonstrated to be one of the state-of-the-art algorithms [26]. In addition, from the results of the OB algorithm, the set of all equilibrated routes can be easily extracted, which will be a precondition for the restriction sensitivity analysis approach in next step of SAB algorithm.The OB algorithm uses theorigin-based approach proportions, αp={αap}, as its main solution variables, where αap represents the proportion of flow that comes from origin p through link a. The approach proportions are updated by shifting flows within the restricting subnetwork, Ap. The details of the OB algorithm can be referred to in Bar-Gera [26]. Using the approach proportions, the OB algorithm is able to store the route flow information with significantly less memory than the route-based algorithms and achieve high-accuracy solutions compared to the link-based algorithms.In the process of the OB algorithm, it starts with trees of minimum cost routes from origins as restricting networks. Then steps of updating restricting network and approach proportions are implemented for each origin separately. An “inner loop” is performed subsequently to accelerate convergence, in which the origin-based link flows are updated while keeping the restricting subnetworks fixed. The inner loop is useful, because the restricting subnetworks tend to stabilize fairly quickly but updating the restricting subnetworks requires more computational effort [26, 27].The OB algorithm for the combined trip distribution and assignment problems requires an additional step to update the O-D flows while keeping the route proportions fixed. Taking Evans’ algorithm as a reference, we modified the original OB algorithm for the solution of the ETDA-VDC model, in which the step size of the O-D flow update is obtained by solving a one-dimension search problem. The proposed algorithm is summarized as shown in Algorithm1.Algorithm 1:
Initialization:
Setn∶=0, determine the initial link cost t0 and destination cost M0 by setting x0=0 and d0=0, respectively.
For each originp∈Z do
(i) Find the tree of minimum cost routes rooting fromp. Let Ap be the minimum route tree. Denote
the minimum route cost to destinationq∈Zp by τpq0, and choose a minimum cost route from pto q.
(ii) Compute the initial variable O-D demands by
dpq0=Opexp⁡[-θ(τpq0+Mq0)]∑kexp⁡[-θ(τpk0+Mk0)],∀q∈Zp
(iii) Forq∈Zp, assign the entire O-D demand epq and dpq0 to the minimum cost route r from p to q, and obtain initial
link flowxap0.
(iv) Update the link costs using the initial link flows.
(v) Initialize the origin-based approach proportionsαap0.
Main Loop:
Given the current variable O-D demanddn obtained in (n-1)th-iteration:
Forn=1 to number of main iterations (IMain)
for eachp in Z do
Update restricting networkAp
Update origin-based approach proportionsαp
end for
Inner Loop:
form=1 to number of inner iterations (Iinner)
for eachp in Z do
Update origin-based approach proportionsαp
end for
end for
UpdateO-Dflows, retain origin-based approach proportions:
Given the origin-based approach proportionsαpn, link flows vn and link costs tn obtained in the steps
above:
(i) For eachq∈Z, compute the destination cost Mqn by O-D demand dpqn and epq.
(ii) Find the set of auxiliary trip demandsd^pqn by solving the following logit distribution model:
d^pqn=Opexp⁡[-θ(τpqn+Mqn)]∑kexp⁡[-θ(τpkn+Mkn)]∀p∈Z,q∈Zp
(iii) Calculate the auxiliary traffic flowv^an on each link a with the approach proportions αapn
(iv) Let(vλ(n),dλ(n))=(1-λ)(vn,dn)+λ(v^n,d^n) and solve the one-dimensional search problem defined
as follows to obtain the step sizeλ*∈[0,1]
min0≤λ≤1⁡F(λ)=∑a∈A∫0vaλ(n)ta(x)dx+(1/θ)∑p∈Z∑q∈Zpdpqλ(n)(ln⁡dpqλ(n)-1)+∑q∈Z∫0Dqλ(n)Mq(y)dy
(v) Set(vn+1,dn+1)∶=(vλ(n),dλ(n)) and check for convergence. Terminate if the convergence
criterion is satisfied; otherwise, update total link flows and link costs, setn∶=n+1 and start
a new iteration of the main loop.
End forThe performance of Algorithm1 was testified on different test networks. The characteristics of these networks are obtained from http://www.bgu.ac.il/~bargera/tntp/. The converging processes on two well-known networks are demonstrated in Figure 1, where the relative gap is defined as RG=|Fn-Fn-1|/Fn (Fn is the value of objective function (7) at nth main loop). It shows that the converging of the improved OB algorithm (IOBA) is faster and more stable than Evans’ algorithm (EA) or the original OB algorithm (OBA) steps in Bar-Gera and Boyce [25].Figure 1
Relative gap versus CPU time.
(a)
(b)


## 3.2. Sensitivity Analysis of the ETDA-VDC Model

On the basis of an equilibrium solution of the lower-level model, the derivatives of the lower-level decision variables with respect to the upper-level ones should be produced for linearly approximating the whole bilevel model in next step. The derivative expressions can be obtained by employing the restriction sensitivity analysis approach for the ETDA-VDC model [18]. The restriction approach reduces the original network to a restricting one in which the nonuniqueness difficulty will be overcome, and thus the full derivatives of the link flows in the primal problem can be solved. The necessary results from the sensitivity analysis of the ETDA-VDC model are presented in this section without proof.In the restricting problem the derivatives of the model solutions with respect to the input parameters are derived from(15)∇Ox(ε)=Jx-1-JO,where Jx is an invertible Jacobian matrix with respect to the solution variables x=(fBT,hBT,dT,λT,μT,uT)T, in which fB and hB are the basic variables of the additional and existing route flows, respectively. JO is the Jacobian matrix with respect to the additional zonal productions, O. Here, λ,μ,u are the Lagrange multipliers associated with the conservation constraints (8)–(10). O is referred to as an input parameter for the ETDA-VDC model. It is in terms of the perturbation ε. The expressions of Jx and JO are as follows:(16)Jx=∇f2L∇f,hLOOO-ΛfT∇h,fL∇h2LOO-ΛhTOOO∇d2L-ΦTOIOO-ΦOOOO-ΛhOOOO-ΛfOIOOO,(17)JO=OOOIOOT,where ∇f2L=ΔfT∇vt(v,ε)Δh, ∇h,fL=ΔhT∇vt(v,ε)Δf, ∇h2L=ΔhT∇vt(v,ε)Δh, ∇f,hL=ΔfT∇vt(v,ε)Δh, and ∇d2L=(1/θ)diag(1/dpq)+ΨT∇dM; diag(1/dpq) is a diagonal matrix with 1/dpq as its diagonal elements. The superscript “T” represents the transposed matrix. “I” indicates the identity matrix. The link cost function, t(·), should be strongly monotonically defined in v, which guarantees the uniqueness of the equilibrium solution and is also a precondition for the sensitivity analysis approach. Utilizing the restriction approach, the derivative, ∇Ox(ε), is first conducted in the restricting problem and then ∇Ov, which is defined in the original one, can be obtained by(18)∇Ov=Δf∇OfB+Δh∇OhB.Equation (18) indicates that the variations on the link flows can be just represented by the changes on the basic route flows, which is the rationale of the restriction sensitivity analysis approach.From (15), the inverse of the matrix Jx is required to be worked out. However, the calculation of the inverse matrix is costly for both computation time and storage space in numerical computation for large-scale problems. The dimension of Jx (from (16)) could be about five times of the number of O-D pairs. So the elements of Jx could totally take approximately O(25n4) space, where n is the amount of the traffic zones, and thus even a medium network with hundreds of origins will conduct a Jacobian matrix with over one hundred thousand dimensions. For example, in a small network with about 30 traffic zones (900 O-D pairs), nearly 400 MB RAM should be expended on the storage of the Jacobian matrix (if each element uses a double type). Consequently, the conventional sensitivity analysis approach for the combined model is very difficult to be utilized in the large-scale problems. Besides, the direct calculation for the inverse of a large matrix is pretty inefficient and inaccurate in practice.To avoid the defects of deriving the inverse matrix directly, the usual way to solve the following linear equation as an alternative. Equation (15) which produces the sensitivity results can be rewritten as(19)Jx·∇Ox=-JO.Then, ∇Ox is derived by solving the following series of linear equations:(20)bi=Jx·zi,where bi is the ith column vector in matrix of [-JO] and zi is the ith column vector in ∇Ox. Let ei be a unit vector with one in the ith position and zeros elsewhere. The length of ei is equal to the column number of JO. Thus, (21)bi=000-ei00,-JO=…,bi,…,where 0 in bold is zero vectors associate with the zero block in JO.To solve (20), Jx is premultiplied by the matrices K1 and K2 in sequence and postmultiplied by K1T and K2T simultaneously, which is equivalent to making the elementary transformation of matrix Jx. Consider(22)K1=IIIIΦII,K2=IIIII-∇d2L-1I.Thus,(23)K2K1JxK1TK2T=∇f2L∇f,hLO-ΛfTΦTO-ΛfT∇h,fL∇h2LOO-ΛhTOOO∇d2LOOO-ΦΛfOOOOOO-ΛhOOOO-ΛfOOOO-∇d2L-1.By using the above operations on (20), we obtain(24)K2K1JxK1TK2TK2T-1K1T-1zi=K2K1bi.Obviously, (25)y=K2T-1K1T-1zi,K2K1bi=000-ei00.Therefore, it is equivalent to solving the equation in terms ofy and then deriving zi by(26)zi=K1TK2Ty=y1y2y3-∇d2L-1·y6y4y5y6+ΦT·y4.Since it is only concerned with the first three subvectors ofzi in the sensitivity result, we just need to find the value of y1, y2, y3, and y6. From (24), the following group of equations can be obtained:(27)∇f2L·y1+∇f,hL·y2-ΛfTΦT·y4-ΛfT·y6=0,(28)∇h,fL·y1+∇h2L·y2-ΛhT·y5=0,(29)∇d2L·y3=0,(30)-ΦΛf·y1=-ei,(31)-Λh·y2=0,(32)-Λf·y1-∇d2L-1·y6=0.The above equation systems can be further simplified depending on whether the incidence matrix Λh is square.(1)
IfΛh is square, consider the following.Thus,Λh will be constructed as an identity matrix; the value of y2 and y3 can be computed from (29) and (31), and thereby y2=y3=0. From (32),(33)y6=-∇d2L·Λf·y1.Substituting (33) into (27), we obtain(34)∇f2L+ΛfT∇d2LΛfy1-ΛfTΦT·y4=0.Combined with (30) we get(35)∇f2L+ΛfT∇d2LΛf-ΛfTΦT-ΦΛfOy1y4=0-ei,where(36)∇f2L+ΛfT∇d2LΛf=ΔfTΛfT∇vt(·)∇d2LΔfΛf.Because the matrix ∇vt(·)∇d2L is positive definite, the columns of ΔfΛf are leaner independent from the restriction approach. It can be easily conducted that matrix (36) is invertible. Let(37)∇f2L+ΛfT∇d2LΛf-ΛfTΦT-ΦΛfO=ABTBO.Solving (35), we can get(38)y1=A-1BTBA-1BT-1·ei.To this extent, the only effort to produce the desired derivative results is to calculatey1, in which the inverse of matrices A and BA-1BT should be derived. The dimension of matrix A is equal to the numbers of the basic routes used by the additional travel demand, which is a little more than the number of O-D pairs from our observation in computational experiments. The dimension of matrix BA-1BT is the same as the number of origins and is not large.Consequently, by utilizing (33) and y2=y3=0, the first three subvectors of zi can be obtained, which are separately corresponding to the derivative results, ∇OfB, ∇OhB, and ∇Od, in ∇Ox. The expressions are(39)∇OfB=A-1BTBA-1BT-1,∇OhB=0,∇Od=Λf·∇OfB.Thus, (40)∇Ov=Δf∇OfB.(2)
IfΛh is not square, consider the following.In this case,y2=∇OhB will not be equal to zero. For further simplification, we rewrite Λh as(41)Λh=Λh′I.Let R′ denote the set of routes associated with Λh′ and the flows on these routes are represented by h′. Repeating the derivation from (33) to (35), we can get(42)∇f2L+ΛfT∇d2LΛf∇f,hL′-ΛfTΦTO∇h,fL′∇h2L′OΛh′T-ΦΛfOOOOΛh′TOO×y1y2′y4y5′=00-ei0,where the superscript “′” indicates the matrices and variables corresponding to the routes in R′. Considering that the set of R′ is fairly small for network capacity problems, the equation system is slightly larger than (35), which is also easy to be solved. The derivative of the existing route flows is computed by y2′. Consider(43)∇OfB=y1,(44)∇OhB=y2′,(45)∇Od=Λf·∇OfB.And thus, (46)∇Ov=Δfy1+Δh′y2′.From the above result, it can be noted that the primal problem to calculate the inverse ofthe entire matrixJx is replaced by solving a small system of linear equations.


## 3.3. Restriction of the Equilibrium Network

When the restriction sensitivity analysis approach was proposed [10], the restricting problem was conducted by finding an extreme point of the feasible region of the equilibrium route flow. An equivalent linear programming was provided for solution. However, this method has been verified usually failed, since the linear programming may find a degenerated solution which causes the fact that the reduced restricting problem cannot reserve the essential information and conducts a wrong result [16]. A rectified method is to find a maximum set of the linear independent columns in the coefficient matrix of the following equation system [18]:(47)Δf0Δh0Λf0OOΛh0f0h0=v*d*e.The above equation system defines the feasible region of the equilbirum route flow under unique euqilibrium solution of v* and d*. Note that the notations are distinguishing from the ones in (16). f0and h0 correspond to the routes that carry positive flows, which are the equilibrated (or minimum-cost) routes according to Assumption 2 in Du et al. [18], and can be obtained from the results of the OB algorithm. The superscript “0” is associated with the variables in a reduced problem, which only consists of the equilibrated routes.From the coefficient matrix of (47), the maximum set of the linear independent columns corresponds to the set of the equilibrated and linear independent routes (ELI, denoted by RB) and can be found using the row echelon form of the coefficient matrix. We can utilize the blocks Λf and Λh which are incidence matrices with one in each column (because one route can only serve one O-D). Since there may be more than one equilibrated route between an O-D pair, Λf and Λh are generally not column full rank. Thus, the coefficient matrix can be rewritten as (48)Λ-fOΛ~fOOΛ-hOΛ~hΔ-fΔ-hΔ~fΔ~h=IΛ~fOIOΛ~hΔ-fΔ-hΔ~fΔ~h⟶IOOOΔ~f-Δ-fΛ~fΔ~h-Δ-hΛ~h,where Λ-f and Λ-h are square matrices, which consist of exactly one equilibrated route for each O-D, so they are equal to identity matrices. Thus, the question is reduced to find and eliminate the linear dependent columns in a submatrix. Consider(49)Δ~f-Δ-fΛ~fΔ~h-Δ-hΛ~h.The columns of the above matrix are only related to the alternative routes between every O-D. Therefore, the set RB will include the routes corresponding to Λ-f, Λ-h, and the maximum set of the independent columns in (49).


## 3.4. Solution of the Maximum Trip Production Problem

When the sensitivity results are derived, the SAB algorithm will use this information to represent the implicit relationships betweenv, d, and O. Thus, the bilevel problem is first-order approximated at the given point, O*=f-1(v*,d*). Let v(O*) and d(O*), respectively, denote the solutions to the lower-level model at O*. The relationship can be represented by using the Taylor expansion:(50)vO≈vO*+∇Ov·O-O*,dO≈dO*+∇Od·O-O*,where the derivatives ∇Ov and ∇Od are obtained from the sensitivity analysis of the ETDA-VDC model. Therefore, the bilevel problem can be reformulated as(51)Max⁡∑p∈ZOps.t.∇Ov·O≤C-v*+∇Ov·O*Φ·∇Od·O≤Omax⁡-O--O*+Φ·∇Od·O*Ψ·∇Od·O≤Dmax⁡-D--D*+Ψ·∇Od·O*.The solution of this linear programming can be easily derived using the simplex method. However, because the linear programming problem is just locally approximated, the new solution might be infeasible to the original problem. At an infeasible point of O* the lower-level trip distribution and assignment results may not satisfy the upper-level constraints. In extreme cases, for example, some link flow va* may be much greater than its capacity Ca, which could cause the capacity constraints to fail for any nonnegative O. In consequence, the new linear approximation will have no solution.In consideration of this flaw, we modified the conventional SAB algorithm by restraining the solution to be always located within the feasible region. A trip distribution and assignment step is implemented at the solution of the above linear programming,O*, and then the results are checked with the capacity conditions. If any capacity constraint is violated, the new solution should be updated by a convex combination of the solution from last iteration, O(n), and O*(n). Therefore, let O^ be the solution to the approximate linear problem above; the maximum step size λ=2-k(k=0,1,2,…) is chosen which ensures the capacity constraint is satisfied at any link or traffic zone. The new solution is updated by(52)O(n+1)∶=O(n)+2-λO^n-On.


## 3.5. A Heuristic Solution Process

In summary, the modified SAB algorithm involves an iterative process between the upper-level and lower-level problems and can be summarized as follows.Step
0 (initialization). Determine an initial value of trip production pattern O(0). Set n∶=0.Step  1 (solving lower-level problem). Solve the ETDA-VDC model in lower-level for the given O(n), by which the equilibrium link flows v(n) and O-D demand d(n) are obtained. The set of all equilibrated routes can be derived at the equilibrium point.Step
2 (sensitivity analysis). Calculate the partial derivatives, ∇Ov and ∇Od, using the sensitivity method for the ETDA-VDC model.Step
3 (local linear approximation). Formulate local linear approximations of the upper-level capacity constraints using the derivative information and solve the approximate linear programming problem to produce an auxiliary trip production O^(n).Step
4 (updating solution). Let O(n+1)∶=O(n)+2-λ(O^(n)-O(n)), where λ=0,1,2,…, until the capacity constraints are satisfied by solving the ETDA-VDC problem with O(n+1). Set n:=n+1.Step  5 (convergence criterion). If |Op(n+1)-Op(n)|≤κ for all p∈Z then stop, where κ is a predetermined tolerance. Otherwise, return to Step 1.The above process is a modification of the conventional SAB approach. The flowchart is shown in Figure2. It should be noted that without the step-size adaption in Step 4 the SAB method only works in very small networks but could fail in larger examples. The set of all equilibrated routes is used to apply the restriction approach for sensitivity analysis. In this study, the modified OB algorithm for the ETDA-VDC problem is used for this purpose. In addition, the derivative results of the sensitivity analysis can give a precise local linear approximation to the upper-level capacity constraints, which is very important to make the heuristic search converge to a considerable good solution.Figure 2
Flowchart of sensitivity analysis based algorithm for network capacity model.In addition, given the nonconvexity of the bilevel problem, the SAB algorithm will converge to a local optimal point [2]. However, as shown in Yang et al. [28], for MPEC models, if the upper-level objective function is a linear function of its decision variables, the heuristic algorithm can at least find a noninferior optimal solution. So the SAB algorithm is supposed to be able to find a satisfying solution for the road network capacity problem in practice.


## 4. Numerical Experiments

Numerical experiments are conducted in this section. In the numerical experiments the link cost function employs the Bureau of Public Roads (BPR) function:(53)tava=tafree1+0.15vaCa4.The destination cost function is defined, referring to Chen and Kasikitwiwat [8], as follows:(54)MqDq=kq∑p∈Zqepq+dpqωq-mq,where kq is a scaling factor between demand and service cost; ωq is a dimensionless parameter related to the severity of congestion; and mq represents a fixed attraction at destination q.The experiments were implemented on Intel Core i5 CPU 3.20 GHz, 4 GB RAM, using the Microsoft Windows 7 operating system. All of the coding was carried out in Visual C#. The solution accuracy is measured in each iterationn by the relative error of the trip productions On (RE) as follows:(55)REn=O(n)-O(n-1)O(n).


### 4.1. Convergent Rate

To testify the efficiency of the SAB method, the Sioux Falls and Anaheim networks are used for the solution of the network capacity model. The IEA algorithm [1] is also employed as a reference in the test. The Sioux Falls is an aggregated network with 24 zones, 24 nodes, and 76 links. The Anaheim is a medium network that consists of 38 zones, 416 nodes, and 914 links (http://www.bgu.ac.il/~bargera/tntp/). The exiting O-D matrix is adapted to make the current traffic flow be unsaturated. The parameters of the destination cost function are given by default. The two algorithms, SAB and IEA, are implemented on these two networks to compare their convergence performances, and the results are listed in Table 1. The stopping tolerance is set to be 10−7.Table 1
Computational result test networks.
Test networks
Observations
SAB algorithm
IEA algorithm
Sioux Falls
CPU time (sec)
320
162
Iterations
18
37
F
1
158498.9
126093.8
Anaheim
CPU time (sec)
650
>6502
Iterations
20
>1632
F
1
137853.3
122039.02
1
F is the value of the objective function when the algorithm terminates.
2IEA does not meet the convergence criteria on the Anaheim network.From the tests, although SAB algorithm shows a slower speed to reach a local optimum, it always produces much better solutions compared to IEA algorithm. In the tests on the Anaheim network, IEA has not converged and tended to zigzag after a few cycles of computation. Since IEA can quickly produce the approximate derivatives at current point, the CPU time taken in each iteration is much less than SAB algorithm. However, IEA could easily stop at a nonoptimal solution or not converge. This indicates that the accuracy of the derivative results has a significant influence on the convergence of the heuristic search for the network capacity model. The converging processes of the two algorithms on the two networks are illustrated in Figure3. Note that the conventional SAB algorithm failed in both networks in our tests.Figure 3
Value of objective function versus CPU time.
(a)
(b)


### 4.2. Convergent Stability

Since the characteristic of the heuristic search, little can be said theoretically about the convergence of the SAB algorithm. In this study, we implement SAB algorithm from different start points to observe the converging process. On both the Sioux Falls and Anaheim networks, SAB algorithm can find a sufficient good solution given any arbitrary initial solution. We select the converging processes from five different start points on Anaheim network and plot the results in Figure4. All the computations were implemented in 20 times of iterations. The final solutions from the five computations are shown in Figure 5.Figure 4
Relative errors of trip production versus iteration from five start points on Anaheim.Figure 5
Solutions of additional trip production from five start points on Anaheim.The experiment result in Figure4 indicates that the solutions by our SAB algorithm can reach a sufficient high precision. All the solutions in Figure 5 are very close, which reflects that the SAB algorithm is comparatively stable for the network capacity problem. Consequently, the proposed SAB algorithm is testified to probably converge to a suboptimal solution from any appropriate feasible start point for the network capacity model. Since the bilevel problem is not convex, the suboptimal solution can be regarded as a satisfying approximation. This solution might be quite close to the global optimum according to the aforementioned characteristics of the network capacity model. Since it only takes a few (less than 30) times of iterations to reach the local optimum, the whole computation of the SAB algorithm can be completed in a considerable CPU time, which is quite superior to the global optimization, for example, genetic algorithm. Thus, the modified SAB method could be a valuable tool for the estimation of the road network capacity in practice.


## 4.1. Convergent Rate

To testify the efficiency of the SAB method, the Sioux Falls and Anaheim networks are used for the solution of the network capacity model. The IEA algorithm [1] is also employed as a reference in the test. The Sioux Falls is an aggregated network with 24 zones, 24 nodes, and 76 links. The Anaheim is a medium network that consists of 38 zones, 416 nodes, and 914 links (http://www.bgu.ac.il/~bargera/tntp/). The exiting O-D matrix is adapted to make the current traffic flow be unsaturated. The parameters of the destination cost function are given by default. The two algorithms, SAB and IEA, are implemented on these two networks to compare their convergence performances, and the results are listed in Table 1. The stopping tolerance is set to be 10−7.Table 1
Computational result test networks.
Test networks
Observations
SAB algorithm
IEA algorithm
Sioux Falls
CPU time (sec)
320
162
Iterations
18
37
F
1
158498.9
126093.8
Anaheim
CPU time (sec)
650
>6502
Iterations
20
>1632
F
1
137853.3
122039.02
1
F is the value of the objective function when the algorithm terminates.
2IEA does not meet the convergence criteria on the Anaheim network.From the tests, although SAB algorithm shows a slower speed to reach a local optimum, it always produces much better solutions compared to IEA algorithm. In the tests on the Anaheim network, IEA has not converged and tended to zigzag after a few cycles of computation. Since IEA can quickly produce the approximate derivatives at current point, the CPU time taken in each iteration is much less than SAB algorithm. However, IEA could easily stop at a nonoptimal solution or not converge. This indicates that the accuracy of the derivative results has a significant influence on the convergence of the heuristic search for the network capacity model. The converging processes of the two algorithms on the two networks are illustrated in Figure3. Note that the conventional SAB algorithm failed in both networks in our tests.Figure 3
Value of objective function versus CPU time.
(a)
(b)


## 4.2. Convergent Stability

Since the characteristic of the heuristic search, little can be said theoretically about the convergence of the SAB algorithm. In this study, we implement SAB algorithm from different start points to observe the converging process. On both the Sioux Falls and Anaheim networks, SAB algorithm can find a sufficient good solution given any arbitrary initial solution. We select the converging processes from five different start points on Anaheim network and plot the results in Figure4. All the computations were implemented in 20 times of iterations. The final solutions from the five computations are shown in Figure 5.Figure 4
Relative errors of trip production versus iteration from five start points on Anaheim.Figure 5
Solutions of additional trip production from five start points on Anaheim.The experiment result in Figure4 indicates that the solutions by our SAB algorithm can reach a sufficient high precision. All the solutions in Figure 5 are very close, which reflects that the SAB algorithm is comparatively stable for the network capacity problem. Consequently, the proposed SAB algorithm is testified to probably converge to a suboptimal solution from any appropriate feasible start point for the network capacity model. Since the bilevel problem is not convex, the suboptimal solution can be regarded as a satisfying approximation. This solution might be quite close to the global optimum according to the aforementioned characteristics of the network capacity model. Since it only takes a few (less than 30) times of iterations to reach the local optimum, the whole computation of the SAB algorithm can be completed in a considerable CPU time, which is quite superior to the global optimization, for example, genetic algorithm. Thus, the modified SAB method could be a valuable tool for the estimation of the road network capacity in practice.


## 5. Conclusions

This paper has presented an effective SAB method for the solution of the transportation network capacity problem. An OB algorithm for ETDA-VDC problem is presented to solve the lower-level problem rapidly and efficiently. The primal problem of sensitivity analysis has been simplified, so the inverse calculation of a full-size matrix has been converted to solving corresponding linear systems of equations. To ensure the heuristic search can proceed properly, the solution from the local linear approximation is rectified to be restrained in its feasible region. Based on these improvements, the SAB method is able to be applied to estimate the capacity of the large-scale real networks. The performance of the proposed SAB method has been demonstrated on two experimental networks.Based on the method in this paper, the network capacity model could be used to provide an ideal travel demand pattern for a given transportation system, by which the traffic resources can be utilized in maximum. From the model results, the traffic engineers can evaluate the design scheme of a transportation network more accurately before its construction. The results can also be used as a reference for the layout of land use in the new developing areas. In addition, the traffic flow pattern at full saturation may reflect the potential bottleneck of the given network. Further applications of the network capacity model could be detected based on the proposed estimation method.



---

*Source: 827094-2015-03-22.xml*