- Email: [email protected]
An efficient computation algorithm for area traffic control problem with link capacity expansions
Applied Mathematics and Computation 188 (2007) 1094–1102 www.elsevier.com/locate/amc
An efficient computation algorithm for area traffic control problem with link capacity expansions Suh-Wen Chiou Department of Information Management, National Dong Hwa University, 1, Sec. 2, Da Hsueh Road, Shou-Feng, Hualien 97401, Taiwan
Abstract An area traffic control problem is considered where the set of link capacity expansions and signal setting variables are simultaneously determined. This problem can be formulated as an optimization problem by taking the user equilibrium traffic assignment as a constraint. In this paper, an efficient computation algorithm based on the projected Quasi-Newton method is presented to effectively solve the area traffic control problem. Numerical calculations on a real data road network are conducted. As it shows, the proposed method combining the locally optimal search and global search heuristic achieved substantially better performance than did traditional approaches in solving the area traffic control problem with expansions of link capacity. 2006 Elsevier Inc. All rights reserved. Keywords: Network design; Area traffic control; Projected Quasi-Newton method
1. Introduction An area traffic control problem with expansion of link capacity is considered where the sum of total travel times and investment cost of link capacity expansion is minimized. The decision variables are the set of link capacity expansion and signal settings. Consider the solutions for network design, Abdulaal and LeBlanc [1] were the first ones who proposed the Hooke–Jeeves algorithm to solve a continuous network design problem and tested that algorithm on a medium-sized realistic network. Suwansirikul et al. [10] proposed a simple heuristic called the Equilibrium Decomposed Optimization (EDO) and performed this heuristic on several example networks. The computing efficiency of EDO results from the decomposition of the network design problem into a set of interacting optimization problems and in each iteration only one user equilibrium needs to be calculated when updating the improvements of all links of the network. Using the results from sensitivity analysis for equilibrium network flows [11], Yang and Yagar [12] proposed a sensitivity analysis based method (SAB) to solve a signalized network design problem. Furthermore, Meng et al. [6] transferred the bilevel program of the network design problem into a single level continuously differentiable optimization problem where an augmented Lagrangian method was proposed and numerical tests were conducted. E-mail address: [email protected] 0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.09.129
S.-W. Chiou / Applied Mathematics and Computation 188 (2007) 1094–1102
1095
A network design problem is a problem of finding the optimal design parameters while taking into account the route choice of users. This problem can be formulated as a special case of the Optimization Program with Equilibrium Constraints (OPEC) problem, which has been extensively discussed by Luo et al. [5], Shimizu et al. [8] and Outrata et al. [7]. In this paper, an OPEC program for an area traffic control problem with equilibrium flows subject to link capacity expansion is investigated, among which the coordinated signal timing plan is defined by the common cycle time, the start and duration of greens. The performance index is defined as the sum of a weighted linear combination of rate of delay and number of stops per unit time for all traffic streams and the investment for link capacity expansions. A projected Quasi-Newton method is proposed to effectively explore the local optima over the search space. Numerical calculations on a real data road network are conducted. As it shows, the proposed algorithm combines the locally optimal search and global search heuristic achieved substantially better performance than did those other approaches when solving the area traffic control problem with expansion of link capacity. The rest of the paper is organized as follows. In the next section, notation and problem formulation are given for the area traffic control problem subject to link capacity expansions. In the Section 3, a projected Quasi-Newton method is proposed and the KKT points are identified where a near optimum can be found effectively. In the Section 4, numerical computations are conducted on a real data Sioux Falls network where promising results demonstrate the performance and robustness of the proposed method. Conclusions and discussions for this paper are made in the Section 5.
2. Problem formulations The problem for solving an area traffic control with link capacity expansion subject to travelers’ route choice is formulated as an OPEC problem in the following way. First, a generalized traffic assignment with decision variables of link capacity expansions and signal settings is presented. Next, an OPEC program is presented where the KKT system of equations is considered. Notation used throughout this paper is summarized below. 2.1. Notation G(N, L) a directed road network, where N is the set of nodes and L is the set of links y the vector of link capacity expansions u the vector of upper bound for link capacity expansions h the conversion factor from investment cost to travel times G(y) the vector of investment costs W = (f, h, /) the set of signal setting variables, respectively for the reciprocal of cycle time, start and duration of greens Ka the duration of effective green for link a gjm the minimum green for signal group j at junction m cjlm the clearance time between the end of green for group j and the start of green for incompatible group l at junction m Xm(j, l) a collection of numbers 0 and 1 for each pair of incompatible signal groups at junction m; where Xm(j, l) = 0 if the start of green for signal group j proceeds that of l, and Xm(j, l) = 1, otherwise Da the rate of delay on link a Sa the number of stops per unit time on link a T the travel demands for OD pairs f the vector of path flows q the vector of link flows d the link-path incidence matrix D the OD-path incidence matrix c the vector of link flow travel time C the vector of path flow travel time
1096
S.-W. Chiou / Applied Mathematics and Computation 188 (2007) 1094–1102
2.2. User equilibrium traffic assignment A general formulation for a user equilibrium traffic assignment problem can be expressed in terms of a variational inequality. To find values q such that ct ðqÞðz qÞ P 0
ð1Þ
for all z 2 K = {q : q = df, T = Df, f P 0} where superscript t is matrix transpose operator. The equivalent variational inequality in terms of path flow can be given as follows. To find values f such that C t ðf Þðz f Þ P 0
ð2Þ
0
for all z 2 K = {f : T = Df, f P 0}. 2.3. Parametric variational inequality For user equilibrium traffic assignment problem with link capacity expansions y and signal settings W = (f, h, /), the parametric variational inequality for (1) can be expressed as follows: ct ðq; y; WÞðz qðy; WÞÞ P 0
ð3Þ
for all z 2 K(y, W) = {q(y, W) : q(y, W) = df(y, W), T = Df(y, W), f(y, W) P 0}. Similarly, the parametric variational inequality in terms of path flow for (2) can be expressed as follows: C t ðf ; y; WÞðz f ðy; WÞÞ P 0
ð4Þ
0
for all z 2 K (y, W) = {f(y, W) : T = Df(y, W), f(y, W) P 0}. 2.4. Sensitivity analysis Suppose (f, y, W) solves problem (4), the KKT system, H, of equations of parametric variational inequality (4) can be expressed. Let w = (f, p, l) and the associated Lagrange multipliers p and l, we have H ðw; y; WÞ ¼ 0; p P 0;
ð5Þ
i.e. Cðf ; y; WÞ p Dt l ¼ 0; diagðpÞf ðy; WÞ ¼ 0; Df ðy; WÞ T ¼ 0; p P 0:
ð6Þ
Suppose we work on path with positive flow, i.e. f > 0, p = 0, the KKT system of equations (5) can be reexpressed as Cðf ; y; WÞ Dt l ¼ 0; Df ðy; WÞ T ¼ 0:
ð7Þ
The first ordersensitivity results for (7) with respect to link capacity expansions and signal settings are as fol y lows. Let e ¼ and w = (f, l), it implies W re H ðw; y; WÞ þ rw H re w ¼ 0; ð8Þ t re C rC D where rw H ¼ . Therefore the first order partial derivatives of equilibrium and re H ¼ 0 D 0 flow and associated Lagrange multiplier with respect to the decision variables are of the following form: re w ¼ rw H 1 ðw; y; WÞre H ðw; y; WÞ:
ð9Þ
S.-W. Chiou / Applied Mathematics and Computation 188 (2007) 1094–1102
1097
2.5. OPEC problem with expansion of link capacity For an area traffic control network with expansion of link capacity, a mathematical program can be given. Let P be the performance index which can be expressed via function P0 in terms of link capacity expansion y, signal setting variables W = (f, h, /) and network flow q. An OPEC problem can be expressed as Min
P ¼ P 0 ðy; W; qÞ
subject to
0 6 y 6 u; W 2 P and
y;W;q
ð10Þ q 2 Kðy; WÞ;
where the set P defines the constraints of signal settings W, which can be expressed as system of linear inequalities for minimum green, maximum cycle time and capacity constraints, and K(y, W) defines the solution set of equilibrium flows. Since the equilibrium traffic assignment can be expressed as a KKT system of equations as given in (7), the OPEC program in (10) can be regarded as a standard non-linear problem (11) with respect to (y, W, q) and therefore is solved by normal non-linear programming techniques Min
P ¼ P 0 ðy; W; qÞ
subject to
0 6 y 6 u;
y;W;q
AWt 6 B
ð11Þ and
H ðy; W; qÞ ¼ 0:
Let W 1aD and W 1aS be respectively link-specific weighting factors for the rate of delay and the number of stops per unit time. The details of model (11) is expressed as X Min P ¼ P 0 ðy; W; qÞ ¼ Da W 1aD þ S a W 1aS þ hGa ðy a Þ y;W;q
subject to
a2L
fmin 6 f 6 fmax ; gjm f 6 /jm 6 1;
8j; m;
qa 6 sa Ka ; 8a 2 L; hjm þ /jm þ cjlm f 6 hlm þ Xm ðj; lÞ; 0 6 y 6 u and
ð12Þ j 6¼ l; 8j; l; m;
H ðy; W; qÞ ¼ 0:
The first constraint is for the common cycle time and for each junction m the constraints [2–4] are for the green phase, link capacity and clearance time. 3. Solution method for (11) Applying the first order partial derivatives given in Chiou [2] for the performance index in (11), the gradients of (11) can be presented as follows at (y0, W0, q0) ry P ðy 0 ; W0 Þ ¼ ry P 0 ðy 0 ; W0 ; q0 Þ þ rq P 0 ðy 0 ; W0 ; q0 Þry qðy 0 ; W0 Þ;
ð13Þ
rW P ðy 0 ; W0 Þ ¼ rW P 0 ðy 0 ; W0 ; q0 Þ þ rq P 0 ðy 0 ; W0 ; q0 ÞrW qðy 0 ; W0 Þ:
ð14Þ
3.1. A projected Quasi-Newton method Regarding the solution method for (11), a projected Quasi-Newton method is presented in the following. Theorem 1 (Quasi-Newton method of Broyden–Fletcher–Goldfarb–Shanno updates). Let X1 locally solve (11), and Q1 positive definite symmetric matrix. For k = 1, . . . , n, let Xk+1 = Xk + akdk, where d k ¼ Qk rP t ðX k Þ
ð15Þ
is a descent direction at Xk with ak for problem (11). For k = 1, . . . , n 1, Qk is given as follows: Qkþ1 ¼ Qk
Qk uk utk Qk vk vtk þ t utk Qk uk v k uk
ð16Þ
1098
S.-W. Chiou / Applied Mathematics and Computation 188 (2007) 1094–1102
where uk = Xk+1 Xk and vtk ¼ rP ðX kþ1 Þ rP ðX k Þ. If $P(Xk) 5 0 for k = 1, . . . , n, then Q1, . . . , Qn are symmetric and positive definite so that d1, . . . , dn are descent directions. Proof. For k = 1, Q1 is symmetric and positive definite, then$P1(X1)td1 = $P1(X1)tQ1$P1(X1) < 0, it implies d1 is a descent direction. We assume that the result holds true for k 6 n 1 and then show that it holds true for k + 1. Let x be a non-zero vector such that 2
xt Qkþ1 x ¼ xt Qk x
2
ðxt Qk uk Þ ðxt vk Þ þ t : t uk Qk uk v k uk
ð17Þ 1
1
1
Since Qk is a symmetric positive definite matrix, there exists a positive definite matrix Q2k such that Qk ¼ Q2k Q2k . 1 1 Let p ¼ Q2k x, and q ¼ Q2k uk , then xtQkx = ptp, utk Qk uk ¼ qt q, and xtQkuk = ptq. Substituting in (17), we have 2
xt Qkþ1 x ¼ pt p
2
2
2
ðpt qÞ ðxt vk Þ ðpt pÞðqt qÞ ðpt qÞ ðxt vk Þ þ þ ¼ : t qt q v k uk qt q vtk uk
ð18Þ
By Schwartz inequality,(ptp)(qtq) (ptq)2 P 0. Also, because Qk is positive definite, it implies qt q ¼utk Qk uk > 0; 8uk 6¼ 0, since dk is a descent direction, uk 5 0. To show that xtQk+1x > 0, it suffices to show that vtk uk > 0 and vk 5 0. Since t
vtk uk ¼ ðrP 1 ðX kþ1 Þ rP 1 ðX k ÞÞ ðX kþ1 X k Þ ¼ ðrP 1 ðX kþ1 Þ rP 1 ðX k ÞÞt ak d k
ðsincerP 1 ðX kþ1 Þd k ¼ 0Þ
t
¼ ak rP 1 ðX k Þ d k t
¼ ak rP 1 ðX k Þ Qk rP 1 ðX k Þ > 0 (since Qk is positive definite and dk is a descent direction with ak > 0). Suppose vk = 0, it will contradicts vtk uk > 0, thus vk 5 0. Therefore xtQk+1x > 0, "x 5 0. It implies Qk+1 is positive definite and follows that t
t
rP 1 ðX kþ1 Þ d kþ1 ¼ rP 1 ðX kþ1 Þ Qkþ1 rP 1 ðX kþ1 Þ < 0;
8rP 1 ðX kþ1 Þ 6¼ 0:
Therefore dk+1 is a descent direction at Xk+1. This completes the proof.
h
The search direction generated by Quasi-Newton method for an unconstrained nonlinear problem is a descent direction which strictly decreases the objective function value at each iteration provided that the corresponding gradient value is not zero. Consider the problem (11), suppose the first partial derivatives for the objective function can be expressed as (13) and (14), the problem (11) can be rewritten as Min y;W
P ¼ P 1 ðy; WÞ
subject to 0 6 y 6 u
and
AWt 6 B:
ð19Þ
In the following, we apply Quasi-Newton method to a linear constraint set as given in (19) by using the concept of gradient projection method proposed by Rosen. Theorem 2 (A projected Quasi-Newton method, when Gk $P1(yk, Wk) 5 0). Let d yk and d W descent k bethe d yk direction determined by (15) and ak the step length which minimize P1 along dk, where d k ¼ . The dW k projection matrix Gk is of the following form: Gk ¼ I M tk ðM k M tk Þ1 M k ;
ð20Þ
Mk is the gradient of active constraints in (19), where the active constraint gradients are linearly independent and thus Mk has full rank. In (19) a sequence of iterates {yk, Wk} can be generated according to y kþ1 ¼ y k þ ak Gk d yk ; Wkþ1 ¼ Wk þ
ak Gk d W k :
ð21Þ ð22Þ
S.-W. Chiou / Applied Mathematics and Computation 188 (2007) 1094–1102
1099
The search direction sk can be determined in the following form: sk ¼ Gk d k :
ð23Þ
Then the sequence of points {yk, Wk} generated by the projected Quasi-Newton method, P 1 ðy k ; Wk Þ > P 1 ðy kþ1 ; Wkþ1 Þ;
k ¼ 0; 1; 2; . . .
ð24Þ
whenever Gk$P1(yk, Wk) 5 0. Proof. Following the results of Theorem 1, we have rP t1 ðy k ; Wk Þd kþ1 ¼ rP t1 ðy k ; Wk ÞrP 1 ðy k ; Wk Þ < 0;
k ¼ 1; 2; 3; . . .
ð25Þ
Multiply Eq. (25) by projection matrix Qk+1, it becomes rP t1 ðy k ; Wk ÞQkþ1 d kþ1 ¼ rP t1 ðy k ; Wk ÞQkþ1 rP 1 ðy k ; Wk Þ ¼ rP t1 ðy k ; Wk ÞQtkþ1 Qkþ1 rP 1 ðy k ; Wk Þ 2
¼ kQkþ1 rP 1 ðy k ; Wk Þk < 0:
ð26Þ
Thus for sufficiently small b, b > 0, we have P 1 ðy k ; Wk Þ > P 1 ððy k ; Wk Þ þ bskþ1 Þ
ð27Þ
Because by definition ak+1 is the step length which minimize P1 along sk+1 from (yk, Wk), it implies P 1 ðy k ; Wk Þ > P 1 ððy k ; Wk Þ þ bskþ1 Þ P P 1 ððy k ; Wk Þ þ akþ1 skþ1 Þ ¼ P 1 ðy kþ1 ; Wkþ1 Þ which completes this proof.
ð28Þ
h
Theorem 3 (A projected Quasi-Newton method, when Gk$P1(yk, Wk) = 0). Following Theorem 2, when Gk$P1(yk, Wk) = 0, if all the Lagrange multipliers corresponding to the active constraint gradients in (19) are positive or zeros, it implies the current (yk, Wk) is a KKT point. Otherwise choose one negative Lagrange multi^ k of the active constraint gradients by deleting the jth row of M ^ k , which cor~j , and construct a new M plier, say l ~j , and make the projection matrix of the following form: responds to the negative component l ^k ¼ I M ^ t Þ1 M ^ k: ^ t ðM ^ kM G k k
ð29Þ
The search direction then can be determined by (23) and the results of Theorem 2 hold. According to the results of Theorems 2 and 3, the search process for solving the problem (19) either can be terminated at a KKT point or a new search direction can be generated. The following lemma is therefore presented. Corollary 1 (Stopping conditions). If (yk, Wk) is a KKT point for problem P1 then the search process may stop; otherwise a new descent direction at (yk, Wk) can be generated according to Theorems 2,3. 3.2. QNEW for solving (11) Consider the solution for problem (11), an efficient algorithm of combining local and global search is given below. 3.2.1. Local search step Local search step is a locally optimal search for a full optimization with respect to the common cycle time, the start and duration of green at each junction which is conducted as follows: Step 1. Start with (yk, Wk) and set index k = 1. Step 2. Solve a user equilibrium traffic assignment problem with link expansion and signal settings, find the first order derivatives by (13) and (14).
1100
S.-W. Chiou / Applied Mathematics and Computation 188 (2007) 1094–1102
Step 3. Use the projected Quasi-Newton method to determine a search direction by (23). Step 4. If Gk$P1(yk, Wk) 5 0, find a new (yk+1, Wk+1) in (21) and (22) by deciding ak the step length, which minimize P1 along sk and letting k k + 1. Go to Step 2. If Gk$P1(yk, Wk) = 0 and all the Lagrange multipliers corresponding to the active constraint gradients are positive or zeros, (yk, Wk) is the KKT point and stop. Otherwise, follow the results of Theorem 3and find a new projection matrix and go to Step 3.
1
2 candidate link for capacity expansion 100
3
12
4
11
5
6
9
8
7
10
16
18
17
13
signalized junction
14
15
23
22
24
21 Fig. 1. Sioux Falls network.
19
20
S.-W. Chiou / Applied Mathematics and Computation 188 (2007) 1094–1102
1101
3.2.2. Global search step Given the cycle time and duration of green at each junction, a global search is implemented in offsets to obtain a better point in another part of the feasible region as a reset point for local search. After a reset point with significant improvement in performance index is found, conduct the local search in Section 3.2.1 again and locate new signal settings until the difference of the values of the performance index between successive iterations is negligible. 4. Numerical calculations In this section, numerical experiments are conducted on a real data Sioux Falls network with selected signalized junctions as shown in Fig. 1. The proposed QNEW, the non-optimal iterative optimization assignment (IOA) and the sensitivity analysis based (SAB) methods are conducted with two distinct sets of initials for area traffic control problem subject to link capacity expansions. The non-signalized link travel time and link investment cost functions used are adopted from Suwansirikul et al. [9, pp. 261–262] where the convex investment function form is adopted from Abdulaal and LeBlanc [1, p. 28], together with the data input details. Computational results are shown in Table 1 where the results for system optimal solutions (SO for short) are computed on the basis of system optimizations (see Chiou [3]). As shown in Fig. 1 10 sets of signalized junctions and candidate links for capacity expansions are taken into account in the following illustration. Performance index is expressed in terms of vehicles while the conversion value is referred to that used in the work of Suwansirikul et al. As it is observed in Table 1, the multiple local optima exist due to the non-convexity of the area traffic control network and evidently each method leads to a different solution.
Table 1 Computational results for Sioux Falls network Variable/algorithm
Initial value of 1/f (in seconds) Initial value of /jm/f(in seconds) Initial value of ya 1/f /13/f /14/f /15/f /16/f /18/f /19/f /1,10/f /1,11/f /1,12/f /1,16/f y(4,11) y(11,4) y(4,5) y(5,4) y(5,9) y(9,5) y(9,10) y(10,9) y(10,11) y(11,10) PI (in veh) #
1st set
2nd set
IOA
SAB
SO
QNEW
IOA
SAB
SO
QNEW
72 28 6.5 110 43 38 40 39 57 52 49 45 37 47 4.6 3.9 1.2 1.5 2.3 2.4 5.6 4.8 5.4 6.2 99.2 78
72 28 6.5 110 51 50 44 49 50 49 48 51 50 50 3.9 3.8 2.6 2.7 2.6 2.9 4.3 4.4 4.7 4.8 95.34 69
72 28 6.5 110 37 39 46 43 44 51 52 53 39 49 5.4 5.3 1.6 1.7 2.8 2.7 5.7 4.9 4.3 4.3 81.7 50
72 28 6.5 110 55 51 57 52 48 50 51 55 59 52 4.9 5.2 1.9 2.6 2.2 2.8 5.7 5.5 4.7 4.1 82.1 24
118 50 12 132 42 40 41 44 59 55 49 53 44 48 5.7 1.6 5.6 1.6 3.1 3.8 7.6 3.8 7.3 3.6 108.7 80
118 50 12 132 47 49 46 41 56 54 48 54 45 49 3.9 3.8 2.6 2.7 2.9 2.9 4.3 4.4 4.7 4.8 94.4 70
118 50 12 132 52 50 49 48 53 57 51 55 52 51 5.4 5.3 1.8 1.7 2.8 2.7 5.7 4.9 4.8 4.3 81.5 53
118 50 12 132 55 52 57 50 58 55 52 55 56 52 4.9 5.2 1.9 2.1 2.3 2.7 5.7 5.5 4.5 4.7 82 27
Where /jm/f denotes the duration of greens for signal group j at junction m measured in sec and 1/f denotes the common cycle time measured in seconds. PI denotes the performance index value measured in veh-h/h, and # denotes the number of equilibrium traffic assignment solved.
1102
S.-W. Chiou / Applied Mathematics and Computation 188 (2007) 1094–1102
As it is seen in Table 1, the proposed heuristic QNEW achieved the best performance with the value of 82.1 and 82 veh for the two different initials while the non-optimal solutions–IOA, achieved the worst performance with the highest values of 99.2 and 108.7 veh. The relative differences between the values achieved by the proposed QNEW and SO are less than 1% for the two distinct sets of initials while those did by the SAB are over 15% and those did by the IOA are up to 30%. Regarding the computational efforts required in the computation processes for the area traffic control problem subject to link capacity expansion, as it seen from Table 1, the corresponding number of equilibrium assignments solved by the proposed QNEW is far less than other traditional methods. 5. Conclusions and further issues In this paper, we presented an improved heuristic QNEW for an area traffic control network with link capacity expansions. The proposed method QNEW has been illustrated successfully with promising results on the real data Sioux Falls road network as compared to traditional methods by consistently yielding better performance and less computational efforts. Consider further issues like solving the generalized user equilibrium assignments with asymmetric cost mappings are being taken into account and empirically computational experiments are being carried out. Acknowledgement Thanks go to Taiwan National Science Council via grant NSC-95-2416-H-259-014. References [1] M. Abdulaal, L.J. LeBlanc, Continuous equilibrium network design models, Transportation Research 13B (1979) 19–32. [2] S.-W. Chiou, TRANSYT derivatives for area traffic control optimization with network equilibrium flows, Transportation Research 37B (2003) 267–290. [3] S.-W. Chiou, An optimal scheme for toll pricing problem, Applied Mathematics and Computation, in press, doi:10.1016/ j.amc.2006.04.065. [4] T.L. Friesz, R.L. Tobin, H-L. Cho, N.J. Mehta, Sensitivity analysis based algorithms for mathematical programs with variational inequality constraints, Mathematical Programming 48 (1990) 265–284. [5] Z.Q. Luo, J-S. Pang, D. Ralph, Mathematical Programs with Equilibrium Constraints, Cambridge University Press, UK, 1997. [6] Q. Meng, H. Yang, M.G.H. Bell, An equivalent continuously differentiable model and a locally convergent algorithm for the continuous network design problem, Transportation Research Part B 35 (2001) 83–105. [7] J. Outrata, M. Kocvara, J. Zowe, Nonsmooth Approach to Optimization Problems with Equilibrium Constraints: Theory, Applications and Numerical Results, Kluwer Academic Publishers, Boston, 1998. [8] K. Shimizu, Y. Ishizuka, J. Bard, Nondifferentiable and Two-level Mathematical Programming, Kluwer Academic Publishers, Boston, 1997. [9] S. Suh, T.J. Kim, Solving nonlinear bilevel programming models of the equilibrium network design problem: a comparative review, Annals of Operations Research 34 (1992) 203–218. [10] C. Suwansirikul, T.L. Friesz, R.L. Tobin, Equilibrium decomposed optimization: a heuristic for continuous equilibrium network design problem, Transportation Science 21 (1987) 254–263. [11] R.L. Tobin, T.L. Friesz, Sensitivity analysis for equilibrium network flow, Transportation Science 22 (1988) 242–250. [12] H. Yang, S. Yagar, Traffic assignment and signal control in saturated road networks, Transportation Research 29A (1995) 125–139.
An efficient computation algorithm for area traffic control problem with link capacity expansions Suh-Wen Chiou Department of Information Management, National Dong Hwa University, 1, Sec. 2, Da Hsueh Road, Shou-Feng, Hualien 97401, Taiwan
Abstract An area traffic control problem is considered where the set of link capacity expansions and signal setting variables are simultaneously determined. This problem can be formulated as an optimization problem by taking the user equilibrium traffic assignment as a constraint. In this paper, an efficient computation algorithm based on the projected Quasi-Newton method is presented to effectively solve the area traffic control problem. Numerical calculations on a real data road network are conducted. As it shows, the proposed method combining the locally optimal search and global search heuristic achieved substantially better performance than did traditional approaches in solving the area traffic control problem with expansions of link capacity. 2006 Elsevier Inc. All rights reserved. Keywords: Network design; Area traffic control; Projected Quasi-Newton method
1. Introduction An area traffic control problem with expansion of link capacity is considered where the sum of total travel times and investment cost of link capacity expansion is minimized. The decision variables are the set of link capacity expansion and signal settings. Consider the solutions for network design, Abdulaal and LeBlanc [1] were the first ones who proposed the Hooke–Jeeves algorithm to solve a continuous network design problem and tested that algorithm on a medium-sized realistic network. Suwansirikul et al. [10] proposed a simple heuristic called the Equilibrium Decomposed Optimization (EDO) and performed this heuristic on several example networks. The computing efficiency of EDO results from the decomposition of the network design problem into a set of interacting optimization problems and in each iteration only one user equilibrium needs to be calculated when updating the improvements of all links of the network. Using the results from sensitivity analysis for equilibrium network flows [11], Yang and Yagar [12] proposed a sensitivity analysis based method (SAB) to solve a signalized network design problem. Furthermore, Meng et al. [6] transferred the bilevel program of the network design problem into a single level continuously differentiable optimization problem where an augmented Lagrangian method was proposed and numerical tests were conducted. E-mail address: [email protected] 0096-3003/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.09.129
S.-W. Chiou / Applied Mathematics and Computation 188 (2007) 1094–1102
1095
A network design problem is a problem of finding the optimal design parameters while taking into account the route choice of users. This problem can be formulated as a special case of the Optimization Program with Equilibrium Constraints (OPEC) problem, which has been extensively discussed by Luo et al. [5], Shimizu et al. [8] and Outrata et al. [7]. In this paper, an OPEC program for an area traffic control problem with equilibrium flows subject to link capacity expansion is investigated, among which the coordinated signal timing plan is defined by the common cycle time, the start and duration of greens. The performance index is defined as the sum of a weighted linear combination of rate of delay and number of stops per unit time for all traffic streams and the investment for link capacity expansions. A projected Quasi-Newton method is proposed to effectively explore the local optima over the search space. Numerical calculations on a real data road network are conducted. As it shows, the proposed algorithm combines the locally optimal search and global search heuristic achieved substantially better performance than did those other approaches when solving the area traffic control problem with expansion of link capacity. The rest of the paper is organized as follows. In the next section, notation and problem formulation are given for the area traffic control problem subject to link capacity expansions. In the Section 3, a projected Quasi-Newton method is proposed and the KKT points are identified where a near optimum can be found effectively. In the Section 4, numerical computations are conducted on a real data Sioux Falls network where promising results demonstrate the performance and robustness of the proposed method. Conclusions and discussions for this paper are made in the Section 5.
2. Problem formulations The problem for solving an area traffic control with link capacity expansion subject to travelers’ route choice is formulated as an OPEC problem in the following way. First, a generalized traffic assignment with decision variables of link capacity expansions and signal settings is presented. Next, an OPEC program is presented where the KKT system of equations is considered. Notation used throughout this paper is summarized below. 2.1. Notation G(N, L) a directed road network, where N is the set of nodes and L is the set of links y the vector of link capacity expansions u the vector of upper bound for link capacity expansions h the conversion factor from investment cost to travel times G(y) the vector of investment costs W = (f, h, /) the set of signal setting variables, respectively for the reciprocal of cycle time, start and duration of greens Ka the duration of effective green for link a gjm the minimum green for signal group j at junction m cjlm the clearance time between the end of green for group j and the start of green for incompatible group l at junction m Xm(j, l) a collection of numbers 0 and 1 for each pair of incompatible signal groups at junction m; where Xm(j, l) = 0 if the start of green for signal group j proceeds that of l, and Xm(j, l) = 1, otherwise Da the rate of delay on link a Sa the number of stops per unit time on link a T the travel demands for OD pairs f the vector of path flows q the vector of link flows d the link-path incidence matrix D the OD-path incidence matrix c the vector of link flow travel time C the vector of path flow travel time
1096
S.-W. Chiou / Applied Mathematics and Computation 188 (2007) 1094–1102
2.2. User equilibrium traffic assignment A general formulation for a user equilibrium traffic assignment problem can be expressed in terms of a variational inequality. To find values q such that ct ðqÞðz qÞ P 0
ð1Þ
for all z 2 K = {q : q = df, T = Df, f P 0} where superscript t is matrix transpose operator. The equivalent variational inequality in terms of path flow can be given as follows. To find values f such that C t ðf Þðz f Þ P 0
ð2Þ
0
for all z 2 K = {f : T = Df, f P 0}. 2.3. Parametric variational inequality For user equilibrium traffic assignment problem with link capacity expansions y and signal settings W = (f, h, /), the parametric variational inequality for (1) can be expressed as follows: ct ðq; y; WÞðz qðy; WÞÞ P 0
ð3Þ
for all z 2 K(y, W) = {q(y, W) : q(y, W) = df(y, W), T = Df(y, W), f(y, W) P 0}. Similarly, the parametric variational inequality in terms of path flow for (2) can be expressed as follows: C t ðf ; y; WÞðz f ðy; WÞÞ P 0
ð4Þ
0
for all z 2 K (y, W) = {f(y, W) : T = Df(y, W), f(y, W) P 0}. 2.4. Sensitivity analysis Suppose (f, y, W) solves problem (4), the KKT system, H, of equations of parametric variational inequality (4) can be expressed. Let w = (f, p, l) and the associated Lagrange multipliers p and l, we have H ðw; y; WÞ ¼ 0; p P 0;
ð5Þ
i.e. Cðf ; y; WÞ p Dt l ¼ 0; diagðpÞf ðy; WÞ ¼ 0; Df ðy; WÞ T ¼ 0; p P 0:
ð6Þ
Suppose we work on path with positive flow, i.e. f > 0, p = 0, the KKT system of equations (5) can be reexpressed as Cðf ; y; WÞ Dt l ¼ 0; Df ðy; WÞ T ¼ 0:
ð7Þ
The first ordersensitivity results for (7) with respect to link capacity expansions and signal settings are as fol y lows. Let e ¼ and w = (f, l), it implies W re H ðw; y; WÞ þ rw H re w ¼ 0; ð8Þ t re C rC D where rw H ¼ . Therefore the first order partial derivatives of equilibrium and re H ¼ 0 D 0 flow and associated Lagrange multiplier with respect to the decision variables are of the following form: re w ¼ rw H 1 ðw; y; WÞre H ðw; y; WÞ:
ð9Þ
S.-W. Chiou / Applied Mathematics and Computation 188 (2007) 1094–1102
1097
2.5. OPEC problem with expansion of link capacity For an area traffic control network with expansion of link capacity, a mathematical program can be given. Let P be the performance index which can be expressed via function P0 in terms of link capacity expansion y, signal setting variables W = (f, h, /) and network flow q. An OPEC problem can be expressed as Min
P ¼ P 0 ðy; W; qÞ
subject to
0 6 y 6 u; W 2 P and
y;W;q
ð10Þ q 2 Kðy; WÞ;
where the set P defines the constraints of signal settings W, which can be expressed as system of linear inequalities for minimum green, maximum cycle time and capacity constraints, and K(y, W) defines the solution set of equilibrium flows. Since the equilibrium traffic assignment can be expressed as a KKT system of equations as given in (7), the OPEC program in (10) can be regarded as a standard non-linear problem (11) with respect to (y, W, q) and therefore is solved by normal non-linear programming techniques Min
P ¼ P 0 ðy; W; qÞ
subject to
0 6 y 6 u;
y;W;q
AWt 6 B
ð11Þ and
H ðy; W; qÞ ¼ 0:
Let W 1aD and W 1aS be respectively link-specific weighting factors for the rate of delay and the number of stops per unit time. The details of model (11) is expressed as X Min P ¼ P 0 ðy; W; qÞ ¼ Da W 1aD þ S a W 1aS þ hGa ðy a Þ y;W;q
subject to
a2L
fmin 6 f 6 fmax ; gjm f 6 /jm 6 1;
8j; m;
qa 6 sa Ka ; 8a 2 L; hjm þ /jm þ cjlm f 6 hlm þ Xm ðj; lÞ; 0 6 y 6 u and
ð12Þ j 6¼ l; 8j; l; m;
H ðy; W; qÞ ¼ 0:
The first constraint is for the common cycle time and for each junction m the constraints [2–4] are for the green phase, link capacity and clearance time. 3. Solution method for (11) Applying the first order partial derivatives given in Chiou [2] for the performance index in (11), the gradients of (11) can be presented as follows at (y0, W0, q0) ry P ðy 0 ; W0 Þ ¼ ry P 0 ðy 0 ; W0 ; q0 Þ þ rq P 0 ðy 0 ; W0 ; q0 Þry qðy 0 ; W0 Þ;
ð13Þ
rW P ðy 0 ; W0 Þ ¼ rW P 0 ðy 0 ; W0 ; q0 Þ þ rq P 0 ðy 0 ; W0 ; q0 ÞrW qðy 0 ; W0 Þ:
ð14Þ
3.1. A projected Quasi-Newton method Regarding the solution method for (11), a projected Quasi-Newton method is presented in the following. Theorem 1 (Quasi-Newton method of Broyden–Fletcher–Goldfarb–Shanno updates). Let X1 locally solve (11), and Q1 positive definite symmetric matrix. For k = 1, . . . , n, let Xk+1 = Xk + akdk, where d k ¼ Qk rP t ðX k Þ
ð15Þ
is a descent direction at Xk with ak for problem (11). For k = 1, . . . , n 1, Qk is given as follows: Qkþ1 ¼ Qk
Qk uk utk Qk vk vtk þ t utk Qk uk v k uk
ð16Þ
1098
S.-W. Chiou / Applied Mathematics and Computation 188 (2007) 1094–1102
where uk = Xk+1 Xk and vtk ¼ rP ðX kþ1 Þ rP ðX k Þ. If $P(Xk) 5 0 for k = 1, . . . , n, then Q1, . . . , Qn are symmetric and positive definite so that d1, . . . , dn are descent directions. Proof. For k = 1, Q1 is symmetric and positive definite, then$P1(X1)td1 = $P1(X1)tQ1$P1(X1) < 0, it implies d1 is a descent direction. We assume that the result holds true for k 6 n 1 and then show that it holds true for k + 1. Let x be a non-zero vector such that 2
xt Qkþ1 x ¼ xt Qk x
2
ðxt Qk uk Þ ðxt vk Þ þ t : t uk Qk uk v k uk
ð17Þ 1
1
1
Since Qk is a symmetric positive definite matrix, there exists a positive definite matrix Q2k such that Qk ¼ Q2k Q2k . 1 1 Let p ¼ Q2k x, and q ¼ Q2k uk , then xtQkx = ptp, utk Qk uk ¼ qt q, and xtQkuk = ptq. Substituting in (17), we have 2
xt Qkþ1 x ¼ pt p
2
2
2
ðpt qÞ ðxt vk Þ ðpt pÞðqt qÞ ðpt qÞ ðxt vk Þ þ þ ¼ : t qt q v k uk qt q vtk uk
ð18Þ
By Schwartz inequality,(ptp)(qtq) (ptq)2 P 0. Also, because Qk is positive definite, it implies qt q ¼utk Qk uk > 0; 8uk 6¼ 0, since dk is a descent direction, uk 5 0. To show that xtQk+1x > 0, it suffices to show that vtk uk > 0 and vk 5 0. Since t
vtk uk ¼ ðrP 1 ðX kþ1 Þ rP 1 ðX k ÞÞ ðX kþ1 X k Þ ¼ ðrP 1 ðX kþ1 Þ rP 1 ðX k ÞÞt ak d k
ðsincerP 1 ðX kþ1 Þd k ¼ 0Þ
t
¼ ak rP 1 ðX k Þ d k t
¼ ak rP 1 ðX k Þ Qk rP 1 ðX k Þ > 0 (since Qk is positive definite and dk is a descent direction with ak > 0). Suppose vk = 0, it will contradicts vtk uk > 0, thus vk 5 0. Therefore xtQk+1x > 0, "x 5 0. It implies Qk+1 is positive definite and follows that t
t
rP 1 ðX kþ1 Þ d kþ1 ¼ rP 1 ðX kþ1 Þ Qkþ1 rP 1 ðX kþ1 Þ < 0;
8rP 1 ðX kþ1 Þ 6¼ 0:
Therefore dk+1 is a descent direction at Xk+1. This completes the proof.
h
The search direction generated by Quasi-Newton method for an unconstrained nonlinear problem is a descent direction which strictly decreases the objective function value at each iteration provided that the corresponding gradient value is not zero. Consider the problem (11), suppose the first partial derivatives for the objective function can be expressed as (13) and (14), the problem (11) can be rewritten as Min y;W
P ¼ P 1 ðy; WÞ
subject to 0 6 y 6 u
and
AWt 6 B:
ð19Þ
In the following, we apply Quasi-Newton method to a linear constraint set as given in (19) by using the concept of gradient projection method proposed by Rosen. Theorem 2 (A projected Quasi-Newton method, when Gk $P1(yk, Wk) 5 0). Let d yk and d W descent k bethe d yk direction determined by (15) and ak the step length which minimize P1 along dk, where d k ¼ . The dW k projection matrix Gk is of the following form: Gk ¼ I M tk ðM k M tk Þ1 M k ;
ð20Þ
Mk is the gradient of active constraints in (19), where the active constraint gradients are linearly independent and thus Mk has full rank. In (19) a sequence of iterates {yk, Wk} can be generated according to y kþ1 ¼ y k þ ak Gk d yk ; Wkþ1 ¼ Wk þ
ak Gk d W k :
ð21Þ ð22Þ
S.-W. Chiou / Applied Mathematics and Computation 188 (2007) 1094–1102
1099
The search direction sk can be determined in the following form: sk ¼ Gk d k :
ð23Þ
Then the sequence of points {yk, Wk} generated by the projected Quasi-Newton method, P 1 ðy k ; Wk Þ > P 1 ðy kþ1 ; Wkþ1 Þ;
k ¼ 0; 1; 2; . . .
ð24Þ
whenever Gk$P1(yk, Wk) 5 0. Proof. Following the results of Theorem 1, we have rP t1 ðy k ; Wk Þd kþ1 ¼ rP t1 ðy k ; Wk ÞrP 1 ðy k ; Wk Þ < 0;
k ¼ 1; 2; 3; . . .
ð25Þ
Multiply Eq. (25) by projection matrix Qk+1, it becomes rP t1 ðy k ; Wk ÞQkþ1 d kþ1 ¼ rP t1 ðy k ; Wk ÞQkþ1 rP 1 ðy k ; Wk Þ ¼ rP t1 ðy k ; Wk ÞQtkþ1 Qkþ1 rP 1 ðy k ; Wk Þ 2
¼ kQkþ1 rP 1 ðy k ; Wk Þk < 0:
ð26Þ
Thus for sufficiently small b, b > 0, we have P 1 ðy k ; Wk Þ > P 1 ððy k ; Wk Þ þ bskþ1 Þ
ð27Þ
Because by definition ak+1 is the step length which minimize P1 along sk+1 from (yk, Wk), it implies P 1 ðy k ; Wk Þ > P 1 ððy k ; Wk Þ þ bskþ1 Þ P P 1 ððy k ; Wk Þ þ akþ1 skþ1 Þ ¼ P 1 ðy kþ1 ; Wkþ1 Þ which completes this proof.
ð28Þ
h
Theorem 3 (A projected Quasi-Newton method, when Gk$P1(yk, Wk) = 0). Following Theorem 2, when Gk$P1(yk, Wk) = 0, if all the Lagrange multipliers corresponding to the active constraint gradients in (19) are positive or zeros, it implies the current (yk, Wk) is a KKT point. Otherwise choose one negative Lagrange multi^ k of the active constraint gradients by deleting the jth row of M ^ k , which cor~j , and construct a new M plier, say l ~j , and make the projection matrix of the following form: responds to the negative component l ^k ¼ I M ^ t Þ1 M ^ k: ^ t ðM ^ kM G k k
ð29Þ
The search direction then can be determined by (23) and the results of Theorem 2 hold. According to the results of Theorems 2 and 3, the search process for solving the problem (19) either can be terminated at a KKT point or a new search direction can be generated. The following lemma is therefore presented. Corollary 1 (Stopping conditions). If (yk, Wk) is a KKT point for problem P1 then the search process may stop; otherwise a new descent direction at (yk, Wk) can be generated according to Theorems 2,3. 3.2. QNEW for solving (11) Consider the solution for problem (11), an efficient algorithm of combining local and global search is given below. 3.2.1. Local search step Local search step is a locally optimal search for a full optimization with respect to the common cycle time, the start and duration of green at each junction which is conducted as follows: Step 1. Start with (yk, Wk) and set index k = 1. Step 2. Solve a user equilibrium traffic assignment problem with link expansion and signal settings, find the first order derivatives by (13) and (14).
1100
S.-W. Chiou / Applied Mathematics and Computation 188 (2007) 1094–1102
Step 3. Use the projected Quasi-Newton method to determine a search direction by (23). Step 4. If Gk$P1(yk, Wk) 5 0, find a new (yk+1, Wk+1) in (21) and (22) by deciding ak the step length, which minimize P1 along sk and letting k k + 1. Go to Step 2. If Gk$P1(yk, Wk) = 0 and all the Lagrange multipliers corresponding to the active constraint gradients are positive or zeros, (yk, Wk) is the KKT point and stop. Otherwise, follow the results of Theorem 3and find a new projection matrix and go to Step 3.
1
2 candidate link for capacity expansion 100
3
12
4
11
5
6
9
8
7
10
16
18
17
13
signalized junction
14
15
23
22
24
21 Fig. 1. Sioux Falls network.
19
20
S.-W. Chiou / Applied Mathematics and Computation 188 (2007) 1094–1102
1101
3.2.2. Global search step Given the cycle time and duration of green at each junction, a global search is implemented in offsets to obtain a better point in another part of the feasible region as a reset point for local search. After a reset point with significant improvement in performance index is found, conduct the local search in Section 3.2.1 again and locate new signal settings until the difference of the values of the performance index between successive iterations is negligible. 4. Numerical calculations In this section, numerical experiments are conducted on a real data Sioux Falls network with selected signalized junctions as shown in Fig. 1. The proposed QNEW, the non-optimal iterative optimization assignment (IOA) and the sensitivity analysis based (SAB) methods are conducted with two distinct sets of initials for area traffic control problem subject to link capacity expansions. The non-signalized link travel time and link investment cost functions used are adopted from Suwansirikul et al. [9, pp. 261–262] where the convex investment function form is adopted from Abdulaal and LeBlanc [1, p. 28], together with the data input details. Computational results are shown in Table 1 where the results for system optimal solutions (SO for short) are computed on the basis of system optimizations (see Chiou [3]). As shown in Fig. 1 10 sets of signalized junctions and candidate links for capacity expansions are taken into account in the following illustration. Performance index is expressed in terms of vehicles while the conversion value is referred to that used in the work of Suwansirikul et al. As it is observed in Table 1, the multiple local optima exist due to the non-convexity of the area traffic control network and evidently each method leads to a different solution.
Table 1 Computational results for Sioux Falls network Variable/algorithm
Initial value of 1/f (in seconds) Initial value of /jm/f(in seconds) Initial value of ya 1/f /13/f /14/f /15/f /16/f /18/f /19/f /1,10/f /1,11/f /1,12/f /1,16/f y(4,11) y(11,4) y(4,5) y(5,4) y(5,9) y(9,5) y(9,10) y(10,9) y(10,11) y(11,10) PI (in veh) #
1st set
2nd set
IOA
SAB
SO
QNEW
IOA
SAB
SO
QNEW
72 28 6.5 110 43 38 40 39 57 52 49 45 37 47 4.6 3.9 1.2 1.5 2.3 2.4 5.6 4.8 5.4 6.2 99.2 78
72 28 6.5 110 51 50 44 49 50 49 48 51 50 50 3.9 3.8 2.6 2.7 2.6 2.9 4.3 4.4 4.7 4.8 95.34 69
72 28 6.5 110 37 39 46 43 44 51 52 53 39 49 5.4 5.3 1.6 1.7 2.8 2.7 5.7 4.9 4.3 4.3 81.7 50
72 28 6.5 110 55 51 57 52 48 50 51 55 59 52 4.9 5.2 1.9 2.6 2.2 2.8 5.7 5.5 4.7 4.1 82.1 24
118 50 12 132 42 40 41 44 59 55 49 53 44 48 5.7 1.6 5.6 1.6 3.1 3.8 7.6 3.8 7.3 3.6 108.7 80
118 50 12 132 47 49 46 41 56 54 48 54 45 49 3.9 3.8 2.6 2.7 2.9 2.9 4.3 4.4 4.7 4.8 94.4 70
118 50 12 132 52 50 49 48 53 57 51 55 52 51 5.4 5.3 1.8 1.7 2.8 2.7 5.7 4.9 4.8 4.3 81.5 53
118 50 12 132 55 52 57 50 58 55 52 55 56 52 4.9 5.2 1.9 2.1 2.3 2.7 5.7 5.5 4.5 4.7 82 27
Where /jm/f denotes the duration of greens for signal group j at junction m measured in sec and 1/f denotes the common cycle time measured in seconds. PI denotes the performance index value measured in veh-h/h, and # denotes the number of equilibrium traffic assignment solved.
1102
S.-W. Chiou / Applied Mathematics and Computation 188 (2007) 1094–1102
As it is seen in Table 1, the proposed heuristic QNEW achieved the best performance with the value of 82.1 and 82 veh for the two different initials while the non-optimal solutions–IOA, achieved the worst performance with the highest values of 99.2 and 108.7 veh. The relative differences between the values achieved by the proposed QNEW and SO are less than 1% for the two distinct sets of initials while those did by the SAB are over 15% and those did by the IOA are up to 30%. Regarding the computational efforts required in the computation processes for the area traffic control problem subject to link capacity expansion, as it seen from Table 1, the corresponding number of equilibrium assignments solved by the proposed QNEW is far less than other traditional methods. 5. Conclusions and further issues In this paper, we presented an improved heuristic QNEW for an area traffic control network with link capacity expansions. The proposed method QNEW has been illustrated successfully with promising results on the real data Sioux Falls road network as compared to traditional methods by consistently yielding better performance and less computational efforts. Consider further issues like solving the generalized user equilibrium assignments with asymmetric cost mappings are being taken into account and empirically computational experiments are being carried out. Acknowledgement Thanks go to Taiwan National Science Council via grant NSC-95-2416-H-259-014. References [1] M. Abdulaal, L.J. LeBlanc, Continuous equilibrium network design models, Transportation Research 13B (1979) 19–32. [2] S.-W. Chiou, TRANSYT derivatives for area traffic control optimization with network equilibrium flows, Transportation Research 37B (2003) 267–290. [3] S.-W. Chiou, An optimal scheme for toll pricing problem, Applied Mathematics and Computation, in press, doi:10.1016/ j.amc.2006.04.065. [4] T.L. Friesz, R.L. Tobin, H-L. Cho, N.J. Mehta, Sensitivity analysis based algorithms for mathematical programs with variational inequality constraints, Mathematical Programming 48 (1990) 265–284. [5] Z.Q. Luo, J-S. Pang, D. Ralph, Mathematical Programs with Equilibrium Constraints, Cambridge University Press, UK, 1997. [6] Q. Meng, H. Yang, M.G.H. Bell, An equivalent continuously differentiable model and a locally convergent algorithm for the continuous network design problem, Transportation Research Part B 35 (2001) 83–105. [7] J. Outrata, M. Kocvara, J. Zowe, Nonsmooth Approach to Optimization Problems with Equilibrium Constraints: Theory, Applications and Numerical Results, Kluwer Academic Publishers, Boston, 1998. [8] K. Shimizu, Y. Ishizuka, J. Bard, Nondifferentiable and Two-level Mathematical Programming, Kluwer Academic Publishers, Boston, 1997. [9] S. Suh, T.J. Kim, Solving nonlinear bilevel programming models of the equilibrium network design problem: a comparative review, Annals of Operations Research 34 (1992) 203–218. [10] C. Suwansirikul, T.L. Friesz, R.L. Tobin, Equilibrium decomposed optimization: a heuristic for continuous equilibrium network design problem, Transportation Science 21 (1987) 254–263. [11] R.L. Tobin, T.L. Friesz, Sensitivity analysis for equilibrium network flow, Transportation Science 22 (1988) 242–250. [12] H. Yang, S. Yagar, Traffic assignment and signal control in saturated road networks, Transportation Research 29A (1995) 125–139.