AN MPCC FORMULATION AND ITS SMOOTH SOLUTION ALGORITHM FOR CONTINUOUS NETWORK

Continuous network design problem (CNDP) is searching for a transportation network configuration to minimize the sum of the total system travel time and the investment cost of link capacity expansions by considering that the travellers follow a traditional Wardrop user equilibrium (UE) to choose their routes. In this paper, the CNDP model can be formulated as mathematical programs with complementarity constraints (MPCC) by describing UE as a non-linear complementarity problem (NCP). To address the difficulty resulting from complementarity constraints in MPCC, they are substituted by the Fischer-Burmeister (FB) function, which can be smoothed by the introduction of the smoothing parameter. Therefore, the MPCC can be transformed into a well-behaved non-linear program (NLP) by replacing the complementarity constraints with a smooth equation. Consequently, the solver such as LINDOGLOBAL in GAMS can be used to solve the smooth approximate NLP to obtain the solution to MPCC for modelling CNDP. The numerical experiments on the example from the literature demonstrate that the proposed algorithm is feasible.


INTRODUCTION
The network design problem (NDP) is seeking of a transportation network configuration that minimizes some objective functions, subject to a traditional Wardrop user equilibrium (UE) as the constraints.The continuous NDP (CNDP) has become one of the most computationally intensive problems in the transportation field [1].The CNDP is to determine how to expand the link capacity to minimize the total system travel cost while the users follow the Wardrop's first principle of traffic equilibrium to choose their routes, i.e. no user can decrease their travel time by a unilateral change of route at equilibrium [2].The measurement of the system performance can be described as the sum of total system travel time and the investment cost to expand the link capacity.Due to different objectives for formulating CNDP, it has been modelled as a bi-level programming problem with the upper level (a non-linear programming problem to minimize the total system cost or maximize the social surplus) and the lower level (a UE problem to account for the users' route choice behaviour).Allsop [3] firstly proposed the algorithm for addressing CNDP, and subsequently CNDP has been continuously studied by many researchers during the last five decades.Lots of related publications have grown over time including reviews by Yang and Bell [1] and Farahani et al. [4].Various algorithms have been proposed to solve CNDP (see Table 1 for details).
In this paper, the CNDP's objective is to minimize the sum of the total system travel time and the investment cost by expanding the link capacity, while route choice behaviour of travellers follows UE, which is described by non-linear complementarity problem (NCP).Thus, the CNDP model can be formulated as mathematical programs with complementarity constraints (MPCC).However, solving MPCC is a hard task because the Mangasarian Fromovitz constraint qualification (MFCQ) does not hold at any feasible point [29,30].To circumvent these problems, some algorithmic approaches have focused on avoiding this formulation.Subsequently, researchers have developed special-purpose algorithms for MPCC such as the branchand-bound method [31], the implicit non-smooth approach [32], the piece-wise sequential quadratic programming (SQP) method [33], and the perturbation and penalization approach [34] analysed in Ref. [35].Recently, some exciting new developments have

AN MPCC FORMULATION AND ITS SMOOTH SOLUTION ALGORITHM FOR CONTINUOUS NETWORK DESIGN PROBLEM
programs.Lin and Fukushima [39] have developed a hybrid approach with active set identification to compute a solution or a point with some kind of stationarity by solving a finite number of non-linear programs.For the success of NLP solvers, Leyffer [40] relaxed equivalent condition to replace the usual complementarity condition.Through modelling CNDP as MPCC, Ban et al. [14] relaxed the strict complementarity condition by a relaxation parameter.The relaxed NLP was solved by the existing NLP solvers when this parameter was progressively reduced.Moreover, the relaxation scheme proposed in Ref. [14] can guarantee to solve the original MPCC successfully under certain conditions [41,42].In addition, Fletcher and Leyffer [43] considered solving MPCC as NLP using standard NLP solvers.They demonstrated the numerical experience on a large demonstrated that the gloomy prognosis about the use of transforming MPCC to a well-behaved NLP may have been premature.Several algorithms are proposed for solving MPCC by transforming it to a well-behaved NLP.In particular, Fukushima and Pang [36] considered a smoothing continuation method for the mathematical programming with equilibrium constraints (MPEC).And, under MPEC-linear independence constraint qualification and the asymptotic weak non-degeneracy, they proved that an accumulation point of KKT points is a B-stationary point of the original problem when it satisfies the second-order necessary conditions for the perturbed problems.Subsequently, similar schemes were presented by Scholtes [35] and Lin and Fukushima [37,38].However, these methods need to solve an infinite sequence of non-linear Next, some assumptions used in this paper are presented as follows [6]: 1) The link travel time function t ij (v ij ,y i ),(i,j)!A is strictly increasing and continuously differentiable with respect to the link flow v ij ,(i,j)!A, for any fixed link capacity expansion y ij ,(i,j)!A. 2) The link travel time function t ij (v ij ,y i ),(i,j)!A and , , , ^ĥ h are all continuous with respect to (v ij ,y ij ).
3) The capacity expansion cost function g ij (y ij ),(i,j)!A is continuously differentiable with respect to y ij .

Reformulation of traffic assignment problem
Using the above notations, the set containing all feasible flow distributions for the network, Ω, in terms of link-flow can be described as follows: Ω is a bounded polyhedron because it is comprised of a set of linear equality constraints [21].Here, travellers are assumed to use the Wardrop's UE to describe their route choice behaviour in the transportation network [47].Within each OD pair, a traveller chooses a route to minimize their travel cost.Here, the NCP is used to represent this UE condition.In the link-flow feasible region, Ω, the NCP problem describing UE is to find v and ρ w ("node potential" in Ref. [48]) such that , , , , collection of MPCC test problems, called by MacMPEC, to indicate the suitability of SQP methods for solving MPCC and its out-performance over interior-point solvers regarding speed and reliability.In this paper, the complementarity constraint in MPCC is substituted by a non-smooth equation U(v,x w ,y,ρ w )=0 using Fischer-Burmeister (FB) function { FB proposed in Ref. [44].
Then, the Φ is smoothed by introducing a parameterization W(v,x w ,y,ρ w ,θ) that is differentiable if the scalar θ is non-zero but coincides with Φ when θ =0 [45].
Consequently, the standard NLP solvers such as LIN-DOGLOBAL in GAMS [46] can solve the problem reliably and efficiently.Finally, the numerical experiments on the example from the references demonstrate the feasibility of our model and algorithm.This paper is organized as follows.An MPCC formulation describing CNDP model is presented in Section 2. Section 3 discusses the algorithm for the proposed MPCC model.Section 4 provides numerical experiments on the example from the references to show the feasibility of the MPCC model with the algorithm.Finally, conclusion is given in Section 5.

PROBLEM FORMULATION 2.1 Notations
In this section, the notations are stated: G=(N,A) -transportation network with nodes and links N -set of nodes in G, where N={1,2,...,n} and n denotes the node A -set of links in G, where (i,j)!A denotes the link, i,j!N W -set of OD pairs, and w!W denotes the OD pair y ij -incremental capacity on expanded link (i,j)!A, and y=(y ij ),(i,j)!A denotes the incremental capacity vector l ij , u ij -lower and upper bounds for capacity expansion of link (i,j)!A, respectively.and l=(l ij ), u=(u ij ),(i,j)!A g ij (y ij ) -cost of incremental capacity, y ij , on expanded link (i,j)!A, and g=(g ij ),(i,j)!A d w -travel demand between the OD pair w!W, and the OD demands are given and fixed in this paper x ij w -flow of link (i,j)!A on the OD pair w!W, and -node-link incidence matrix associated with the network, where where z=(v,x w ,y,ρ w ) Now, let us focus on solving the MPCC Problem 6.

SOLUTION ALGORITHM
Here, function {:R 2 →R is called an NCP-function if {(a,b)=0,ab=0,a≥0,b≥0.One popular choice of an NCP-function is the FB function [44]: ( , ) The FB function has many interesting properties that { FB is a convex NCP-function and differentiable on R 2 expect the origin, and { 2 FB is continuously differentiable on R 2 [50].FB function has been used to solve the traffic equilibrium problem [51,52].By using { FB , the complementarity constraint 0≤r ┴ s≥0 can be equivalently transformed into the following non-smooth equation: , .r s r s r s for FB function { FB , the { FB (r,s) is smoothed by introducing a parameterization W(r,s,θ) that is differentiable if the scalar θ is non-zero but coincides with Φ when θ=0 [45].The introduction of the smoothing parameter θ has three consequences [53]: non-smooth problems are transformed into smooth problems, except when θ=0; well-posedness can be improved in the sense that feasibility and constraint qualifications, hence stability, are often more likely to be satisfied for all values of θ; and the solvability of quadratic approximation problems is improved.Therefore, the smooth and approximate equation for the complementarity constraint is obtained as follows: , , r s r s r s We define Ω(θ)={(r,s)|W(r,s,θ)=0}.We analyse the convergence of the smoothing perturbation-based approach by demonstrating the convergence of Ω(θ) to Ω(0) as θ→0.Theorem: For Ω(θ)={(r,s,)|W(r,s,θ)=0}, we have lim 0 Proof: For any , , limsup r s Thus, we have , , , r s r s r s r s 0 0 where symbol " ┴ " is the "perp" operator such that a ┴ b,a T b.The complementarity constraint [30] requires a product of two non-negative variables to be zero, consequently making their values complementary, i.e. when one variable is positive, the other must be zero.It is also possible for both variables to be zero, a case in which the complementarity condition does not hold strictly [49].

The MPCC formulation for the continuous network design problem
In this research, the objective of the CNDP is to minimize the sum of the total system travel time and the investment cost to expand the link capacity, while route choice behaviour of travellers follows UE described by NCP.The CNDP model can be formulated as the following MPCC: , For simplicity, the following notations are introduced.
Step 4: Solution Report The optimal solution is (v * ,x w ,y * ,ρ w * )=z * and the objective function value is f * =f(z k ).

NUMERICAL EXAMPLES 4.1 The 16-link network example
In this subsection, the 16-link network, which was first presented by Harker and Friesz [54], is chosen as the numerical example.The network has been extensively tested in the CNDP literatures such as Ref. [6,8,12,13,14,20,21].As shown in Figure 1, the 16link network consists of two OD pairs, six nodes and links.All input information by Suwansirikul et al. [6] for the test network is presented in Table 2.The travel demand of (1,6) is assumed to be d, which is half of that for (6,1).
In this paper, the numerical experiments consist of three parts.In the first part, we compare the optimal solutions and the objective function values by our proposed algorithm with those by other algorithms in literature.In the second part, we present the results for the test network with different OD demand to make further test on our model and algorithm.In the last part, we consider the total investment cost for link capacity from θ→0 and (r k ,s k )→ (r,s).That is, (r,s)!Ω(0).Therefore we have limsup 0 For any (r,s)!Ω(0), let I + ={i|r i z>}, J + ={i|s i >0}, I 0 ={1,2,...,|A|}\(I + ,J + ), where |A| is the number of links in the transportation network G.For any θ>0, (r i (θ), s i (θ)) is defined by Then , , r s r s r s 0 Obviously ((r(θ),s(θ))→(r,s) and this implies that .liminf 0 Therefore Ω(θ)→ Ω(0) as θ→0.
Hence, the solution to the problem (7) converges to the solution the model (3) as θ→0.
Consequently, the standard NLP solvers to obtain the global solution such as BARON, COIN-OR, LINDO-GLOBAL, LGO, MSNLP and OQNLP in GAMS [46] can be used to deal with Problem 7 reliably and efficiently.The steps of algorithm to solve Model 3 are as follows: Step 1: initialization The parameters are set as follows: θ 0 >θ,e 1 ,e 2 >0 the iteration limit M, l!(0,1).
Step 2: Major Iteration Solve the current relaxed problem (7) , stop and go to Step 4.
algorithm nearly equals that of RELAX.The difference in the objective function values between our algorithm and RELAX is less than 0.02%.
In the second part of the experiments, the model is solved under different OD demand cases.The link capacity expansion with different d of CNDP is shown in Table 4.The upper bound of capacity expansion for each link under these cases are also different.
It is obvious that the objective function values increase with the increase of d.From Table 4, it can be seen that the links y 6 and y 16 are the most necessary ones, whose link capacity should be enhanced.And, y 6 reaches the largest link capacity expansion when d=25, 30,35,40,45,50.The capacity improvement on link y 16 reaches the largest link capacity expansion when d=10,15,30.
In the third part of the experiment, the objective function only considers the total travel times, i.e. , , while the total investment cost to expand the link capacity is considered as a constraint, i.e.

g y B
, , where B is the value of predetermined budget for total investment on link capacity expansions.Therefore, the model is solved with different values of B. Table 5 presents the link capacity expansions for the different B of CNDP.The upper bound of capacity expansion for each link under these cases is identical.
It is obvious that the total travel times decrease as the value of predetermined budget for the total investment cost of link capacity expansions increases.While links y 6 and y 16 are also the targets to enhance link capacity expansions for decreasing the total travel times.
However, links y 6 and y 16 cannot be selected to enhance the link capacity expansions when the value of predetermined budget for total investment cost of link capacity expansions is large enough, such as when B=350,400,450,500. expansion as a constraint by removing it from the objective function to be a constraint function.Then, the model is solved with different values of predetermined budget for total investment on link capacity expansion.In these experiments, a personal computer with Intel Core i5, 3.20 GHz CPU, 4.00 GB RAM, and Windows XP operating system was used for all numerical tests.In addition, GAMS23.5.2 and our familiar LINDOGLOBAL solver was used for solving the smooth approximate problems [46].In numerical experiments, the initial parameters are chosen as follows: θ 0 =1, e 1 =1e-4, e 2 =1e-6, M=15, l=0.2.
In the first part of the experiments, we compare our proposed algorithm with those in the references.A low travel demand level with d=5 and a high travel demand level with d=10 are considered for the tests, which are the same as those in references.Table 3 presents the results of the link capacity expansions under the low and high travel demand cases, respectively obtained by our proposed algorithm and those in the previous studies.Notice that the upper bound of capacity pansion for each link in these two cases is different.The links with zero capacity expansions under all algorithms are not presented in tables for the purpose of space saving.
For the case with low travel demand, the objective function value 199.625 achieved by our algorithm is higher than that by SA, CG, QNEW and LMILP.Among all the previous models, SA, which is regarded to produce the global optimal solution [54], achieves the lowest objective function value 198.10378 .From the value, we can see that the objective function value obtained by our algorithm is very close to that obtained by the SA method.The difference in the objective function values between our algorithm and SA is less than 8%.For the case with high travel demand, the objective function value 522.723 by our algorithm is only higher than that by RELAX.The objective function value by our

Sioux Falls network example
To further test the performance of our proposed algorithm, we have chosen the Sioux Falls network (shown in Figure 2) as another example.This network consists of 24 nodes, 76 links and 552 OD pairs.A subset containing 10 links (16,17,19,20,25,26,29,39,48,74) was chosen for capacity improvement.The link cost function with the corresponding parameters and the OD demand are presented in Tables 6 and 7 [6].The experimental environment and the parameters are adopted as the same as those in Subsection 4.1.We have computed the Sioux Falls network by our proposed algorithm and compared the obtained solution with the corresponding objective function values with those by other algorithms in literature, which are presented in Table 8.Notice that the upper bound of capacity expansion for each link is 25.0.The links with zero capacity expansions under all algorithms are not presented in tables for space saving.

CONCLUSION
In this paper, MPCC has been modelled to describe the CNDP by describing UE as NCP.We have substituted the complementarity constraint in MPCC by a smoothed equation by introducing parameterization W. Therefore, MPCC is transformed into a well-behaved non-linear program (NLP), which can be solved by LINDOGLOBAL in GAMS because W is differentiable if the scalar θ is non-zero but coincides with U when θ=0.The numerical experiments on the 16-link network show that our proposed algorithm can obtain good solutions under the lower and high-level demand cases compared to other algorithms.We have also demonstrated the results about the link capacity enhancement under different OD demand cases and the value of predetermined budget for total investment on link capacity expansions.The Sioux Falls network is also presented to further demonstrate the feasibility of our proposed algorithm.The lower bound of y is set as 0; and, the upper bound of y is set as 25.f is the objective function values from the Ref. [21].f' is the objective function value from reference.B' is the value of the investment on adding link capacity.

ACKNOWLEDGEMENTS
Model 3  for CNDP formulated as MPCC can be transformed into the following smooth and approximate problem:

5 )
Let us rewrite Problem 5 in the following form: for θ=θ k .The obtained solution is z k .Step 3: Stop Condition If k≥M, stop and go to Step 4. Otherwise, compute

Table 1 -
Some algorithms for solving CNDP

Table 2 -
Parameters in the 16-link network

Table 4 -
The results with different OD demand for test network Note: The lower bound of y is set as 0. And, the upper bound of y is set as 2d with different d.Promet -Traffic&Transportation, Vol. 29, 2017, No. 6, 569-580 575

Table 3 -
Comparison of algorithms under different demand for test network

Table 6 -
The link cost function and parameters in Sioux Falls network

Table 7 -
The OD demand in Sioux Falls network

Table 8 -
The obtained results by different algorithms