UNIVERSITY OF

ILLINOIS LIBRARY

AT URBANA-CHAMPAIQN

BOOKSTACKS

va%

CENTRAL CIRCULATION AND BOOKSTACKS

The person borrowing this material is re- sponsible for its renewal or return before the Latest Date stamped below. You may be charged a minimum fee of $75.00 for each non-returned or lost item.

Theft, mutilation, or defacement of library materials can be causes for student disciplinary action. All materials owned by the University of Illinois Library are the property of the State of Illinois and are protected by Article 16B of Illinois Criminal Law and Procedure.

TO RENEW, CALL (217) 333-8400. University of Illinois Library at Urbana-Champaign

SE? 011999

■6 200

When renewing by phone, write new due date below previous due date. L162

BEBR

FACULTY WORKING PAPER NO. 1463

Solving a Real World Highway Network

Design Problem Using Bilevel Linear Programming

Omar Ben-Ayed Charles E. Blair David E. Bo\ce

Hi*

9f

hit

College of Commerce and Business Administration Bureau of Economic and Business Research University of Illinois. Urbana-Champaign

BEBR

FACULTY WORKING PAPER NO. 1463

College of Commerce and Business Administration

University of Illinois at Urbana- Champaign

June 1988

Solving a Real World Highway Network Design Problem Using Bilevel Linear Programming

Omar Ben-Ayed, Graduate Student Department of Business Administration

Charles E. Blair, Professor Department of Business Administration

David E. Boyce Department of Civil Engineering, University of Illinois

Digitized by the Internet Archive

in 2011 with funding from

University of Illinois Urbana-Champaign

http://www.archive.org/details/solvingrealworld1463bena

SOLVING A REAL WORLD HIGHWAY NETWORK DESIGN PROBLEM USING BILEVEL LINEAR PROGRAMMING

OMAR BEN-AYED and CHARLES E. BLAIR

Department of Business Administration University of Illinois at Urbana-Champaign, IL 61820

DAVID E. BOYCE

Department of Civil Engineering University of Illinois at Urbana-Champaign, IL 61801

(May 1988)

This paper is concerned with the solution of a real world Bilevel Linear Program (BLP) of the highway Network Design Problem (NDP), based on empirical data for Tunisia. The problem includes 2,683 variables (2,571 at the lower level) and 820 constraints; it is far beyond the capacity of any existing algorithm. A new algorithm is devised and the solution is found successfully despite the NP-hardness of BLP and the large scale of the problem. The computa- tional difficulties of BLP should not present a barrier against the effective use of such a powerful optimization technique.

Bilevel Linear Programming (BLP) is a mathematical model of an organiza- tional hierarchy in which two decision makers have to choose their strategies from the polyhedron {(x,y): Ax + By < b; x, y > 0). The upper decision maker, who has control over x, makes his decision first, hence fixing x before the lower decision maker selects y [Bialas and Karwan 1982 and 1984, Bard 1983, and, Ben-Ayed 1988] .

The ability of BLP to represent decentralized decision-making problems motivated the development of several applications of the model, including resource allocation [Cassidy 1971, Fortuny-Amat and McCarl 1981], military planning [Bracken and McGill 1973, 1974a and 1974b, Shere and Wingate 1976, Falk 1977] and government policy [Candler and Norton 1977, DeSilva 1978, Bialas et al. 1980, Falk and McCormick 1982, Candler and Townsley 1982 and 1985].

2 Recently, attention has been focused on applying bilevel formulations to a transportation decision making problem referred to as the Network Design Problem (NDP) [LeBlanc and Boyce 1986, Marcotte 1986, Ben-Ayed, Boyce and Blair 1988]. NDP is concerned with the optimal improvement of a transportation network in order to minimize the system total costs.

Most BLP applications are limited to the theoretical formulation, thereby lacking the practical significance of solving real problems. Actually, the overwhelming majority of real-world problems are formulated and solved as single-level (single-objective) programs even when they are virtually bilevel. This is mainly due to the inefficiency associated with BLP algorithms [Ben-Ayed and Blair 1988], which presents a real barrier against the effective use of the model. In trying to motivate further applications of BLP going beyond the theoretical formulations and the small illustrative examples to interesting real world problems, Ben-Ayed, Blair and Boyce [1988] presented a BLP formula- tion of the interregional highway NDP based on empirical data for Tunisia. This paper gives the solution to the problem formulated there. This problem includes 2,683 variables (2,571 of which are lower) and 820 constraints, which makes it too large to be solved by any of the existing algorithms.

Section 1 presents the problem. In Section 2, we develop a new algorithm that takes advantage of the special structure of the problem. There are two reasons for having a BLP formulation; first the user-optimized flow requirement (user-equilibrium), and second the nonconvex improvement functions. The algorithm deals with each of the two lower problems separately; at each iteration we try to find a better compromise with the user, while including the smallest possible number of nonconvex improvement functions to get the exact solution with the minimum computation effort. In Section 3, the solution is presented and the results are analyzed. We conclude in the last section with some suggestions.

1. The Problem to be Solved

In this section we give a short description of the formulation resulting from the empirical study conducted by Ben-Ayed, Blair and Boyce [1988] on the Tunisian interregional highway transportation network. We are concerned with a system of roads connecting the regions of Tunisia. Each region is represented by a centroid-node which is usually the largest city in the region. The traffic entering and leaving a centroid-node is assumed to be that of the whole region; every centroid-node is simultaneously an origin-node and a destination-node. The other nodes in the network are intersections of roads generally correspond- ing to intermediate cities. A road connecting any two cities is called a link and a sequence of links is called a route. We have empirical data on the number of people who want to travel from each region to each other region. Typically, a traveller has a choice of routes and makes his decision based on road conditions. Improvements in various links in the network sill cause some travellers to change their choice of routes. We will not provide details about the empirical data, the construction of the objective functions and the constraints in the formulation; however, a brief discussion on the inclusion of the bilevel presentation is relevant to the object of this paper.

The principal variables in the formulation are the traffic flow and the added capacity resulting from the improvement, both measured in Passenger Car Units per hour (PCU). The designers of a road network wish to find the best combination of improvements to make. Their objective takes into account the cost of improvements and also factors which depend on the pattern of road use, including system travel time, system operating cost and frequency of accidents. The users of the network decide which routes to use primarily based on their individual cost; therefore, their choices do not necessarily coincide with the choices that are optimal for the system [Wardrop 1952]. The system can in- fluence users choices, however, by improving some links and making them more

4 attractive than others. When deciding on these improvements, the system tries to influence the users' preferences in order to minimize total system costs. Thus in our BLP, the upper decision maker chooses the amount of improvement to make on each of the different roads (links in the network) while the lower decision maker chooses how much traffic there will be on each road. A similar approach is used by LeBlanc and Boyce [1986], and Ben-Ayed, Boyce and Blair [1988] .

In addition to the user-optimized flow requirement, BLP formulation allows the analysis of nonconvex piecewise linear functions [Ben-Ayed, Blair and Boyce 1988]. A typical improvement function might be:

f(Z) = 10Z 0 < Z < 200

= 2000 + 3(Z-200), 200 < Z < 500 (1)

= 2900 + 8(Z-500), 500 > 500.

where Z is the amount of improvement of a link and f(Z) is the cost of that improvement. Note that f(Z) is neither convex nor concave, thus making it difficult to analyze. However, the function can be included in the upper objective of a BLP with an appropriate lower objective. The technique is based on the following:

f(Z) = F(Z,W) = 10Z - 7WX + 5W2, where W^MAXfO, Z-200} , W2=MAX( 0 , Z-500] .

The variables Wm are controlled by the lower problem; the lower objective trying to make the Wm as small as possible, and the constraints W^ > Z-200, W2 > Z-500 guarantee that Wm are least upper bounds. This technique can be used for any piecewise linear function.

Two drawbacks should be mentioned. First, the method described does increase the number of lower variables, and this increase may be considerable if there are many piecewise linear functions in the problem. Second, the resulting BLP

5 is typically distinguished by the conflict between its upper and lower objectives; there are many variables which the upper objective wants to be as large as possible, while the lower objective wants them as small as possible. Experience suggests this type of BLP is difficult to solve.

Including user-optimized flow and nonconvex improvement functions, the BLP problem obtained by Ben-Ayed, Blair and Boyce [1988] is synthesized in the following:

MINZ F(Z,W) + G(C)

where [C, C, W, X and Y] solve:

MIN F(W) + G(C)

subject to:

%(¥) = 0

(2)

H2(X,Y) = 0

H3(Z,C,X) > 0

H4(Z,C,X) > 0

H5(Z,W) > 0

H6(Z) > 0

Z, C, C, W, X, Y > 0

where the function F is the improvement cost. The function G is the inter- regional system travel cost. The vector X is the amount of traffic on each link. Any traveller from one origin-node to a destination-node has a choice of different routes to use, corresponding to different sets of links. We use the components of the vector Y to indicate the number of people using each route. These numbers determine the numbers of users of each link, shown by the vector X; the relationship between X and Y is given by the equation constraints H2(X,Y) while the requirement that specified numbers of people go between specified nodes is given by the constraints H^Y) referred to as the conserva-

6 tion of flow conditions. F is the lower objective of nonconvex improvement cost functions. The components of the vector Z are the amounts of improvement on each of the different links. W is the vector of the second-stage variables used in computing the improvement cost functions as described previously in (1). G gives the interregional travel costs on each link for the user. C and C are vectors of second-stage variables giving total travel costs on each link for the system and the user, respectively. H3(Z,C,X) and H^(Z,Z,X) are the ine- quality constraints needed to the formulation of convex travel cost functions and the effect on improvement on decreasing those costs for the system and the user, respectively. Hr(Z,W) correspond to the nonconvex piecewise linear functions constraints. Hg(Z) are the improvement limit constraints. A more detailed description of the variables and the functions is provided in the Appendix A.

In the above formulation (2), the total number of variables (see Appendix A) is 2,683 and the total number of constraints is 820, excluding upper and lower bound constraints (improvement limit constraints and nonnegativity const- raints). Those numbers are partitioned as follows:

Variables Constraints

(excluding bound constraints)

112 (X)

112 (C) 120 (conservation of flow conditions )

112 (C) 112 (link flow definition constraints)

112 (Z) 224 (system travel costs formulation)

104 (W^) 224 (cumulative user costs formulation)

36 (Wo) 68 (convex improvement costs formulation)

2,095 (Y) 72 (nonconvex functions formulation)

2,683 820.

The Z are the only variables controlled by the upper problem; therefore among the 2,683 variables, 2,571 are lower. The nonconvex improvement cost functions all have 3 pieces; two constraints are needed by each link. The convex improvement functions require only one constraint per link because they

all consist of two pieces only.

An ideal solution (lower bound) for any BLP minimization problem can be obtained by ignoring the lower objective function and solving the problem as an LP. If we call (Z°, X°, C°, C°, Y°, W°) the obtained solution, the ideal solution is equal to F(Z ,W ) + G(C ). To get an incumbent solution (upper bound), we substitute Z by Z and solve the lower LP problem. If we call (X , C , C , Y , W ) the obtained solution, the incumbent solution is equal to F(Z ,W ) + G(C ). Let us refer to the gap between incumbent and ideal solutions as:

gap = (incumbent - ideal)/ideal .

The gap for the BLP problem (2) is larger than 111 % [Ben-Ayed 1988]; this high value, mainly due to the nonconvex functions formulation, makes the task of heuristic algorithms extremely difficult, if not impossible. Heuristic al- gorithms, such as the Parametric Complementary Pivot algorithm (PCP) [Bialas- Karwan-Shaw 1980, and Bialas-Karwan 1984] and the Grid Search Algorithm (GSA) [Bard 1983], lead to very unsatisfactory results when applied to BLP problems with combinatorial nature, such as Knapsack problems or nonconvex formulation problems [Ben-Ayed and Blair 1988].

The problem is even harder for branch and bound or related implicit enumera- tion techniques [Falk 1973, Gallo and Ulkiicu 1977, Fortuny-Amat and McCarl

1981, Bialas and Karwan 1982, Papavassilopoulos 1982, Candler and Townsley

1982, Bard and Falk 1982, Bard and Moore 1987]. The Bard-Moore algorithm, which is supposed to be one of the most efficient, does not accept a problem with more than 100 lower variables; our problem has more than 25 times this number.

As for any large NP-hard problem, the use of general algorithms is not efficient. A new algorithm taking advantage of the special structure of the problem has to be devised to solve the problem at hand.

(4)

2. A Proposed Algorithm

Formulation (2) shows that there are two separate lower objective functions. In fact, each of the two objectives has its own constraints that are indepen- dent of the constraints of the other objective, provided that Z are fixed (they are fixed by the upper objective). This can be better seen when the BLP problem (2) is restated as follows:

MINZ>0 F(Z,W) + G(C) (3)

such that {W) solve (4):

MIN F(W)

subject to:

H5(Z,W) > 0

W > 0

and {C, C, X and Y) solve (5):

MIN G(C)

subject to:

H(X,C,C,X,Y) > 0

C, C, X, Y > 0

where H(X,C,C,X,Y) > 0 substitutes the constraints H1(Y) = 0, H2(X,Y) = 0, H3(Z,C,X) > 0, H4(Z,C,X) > 0, and H6(Z) > 0.

The first lower problem (4) corresponds to the nonconvex functions formula- tion and the other lower problem (5) corresponds to the equilibrium problem. Because of the NP-hardness of BLP, it is more convenient to take advantage of the structure of the problem and solve it as two separate smaller problems.

As we mentioned in our discussion of improvement cost functions in the previous section, the BLP corresponding to (4) is difficult. However, the situation is different for the BLP corresponding to the equilibrium problem (5); the travel cost functions of the system and the user are very similar

(5)

9 [Ben-Ayed, Blair and Boyce 1988]. The coefficients of the two functions always have the same sign, and are even fairly proportional; therefore the use of BLP heuristics is safe for this second part.

Our approach to satisfy (5) tries to identify a solution that is good for both the upper and lower objectives by using a convex combination of the two:

tG(C) + (l-t)G(C)

so that we are looking at the problem:

MINZ F(Z,W) + tG(C) + (l-x)G(C)

such that {Z, C, C, X and Y) solve:

H(X,C,C,X,Y) > 0

Z, C, C, X, Y > 0

and such that (W) solve: (6)

MIN F(W)

subject to:

H5(Z,W) > 0

W > 0.

In general, a BLP solution (x ,y ) is called feasible if it belongs to the polyhedron (Ax° + By < b; x , y° > 0}; and satisfies the optimality of the lower solution y to the lower problem when the upper variable x is fixed at x . In our case, for a feasible solution to (6) to be feasible to (2), it must be the case that, for the given choice of the upper variables Z, the values of the lower variables are such as to minimize the lower objective G. This will clearly be the case for t = 0. To make the upper objective as small as pos- sible, we want to find the largest x for which this is the case.

We will assume that if, for some T , the solution to (6) is feasible, then the same will be true for any 0 < t < t . This is not always true, but counter-

10 examples are not easy to construct, and a simulation suggests this is a safe assumption. Granted this assumption, the optimal t, call it x , can be located by a binary search, which is much quicker than a grid search, as used in Bard [1983].

For a given t, (6) is still a BLP with too many lower variables to be solved in a straightforward way. We approach this problem by trying to identify links on which no improvement will be made. For those links, we replace the nonlinear improvement cost function f(Z) by aZ + B, chosen so that aZ + /3 < f(Z) for all Z. We are temporarily assuming that the costs for those links we think are unlikely to be improved are less than they really are. This produces a smaller, more tractable, problem. If the solution to this smaller problem does not use a link with the simplified, cheaper cost it certainly would not want an improve- ment with the real cost. In the case where the small problem does recommend an improvement based on the simplified cost, we must solve a new problem with the real cost as part of the formulation. Thus we solve a sequence of problems in which links with real costs are added when it seems necessary, but not other- wise. Fortunately, in practice the process terminates (no new links added) with a relatively small number of links showing improvement.

Let us call S the set of the links having convex or linear cost improvement cost functions, T the set of the links having nonconvex improvement cost functions, T the subset of T containing all links for which Z are different from 0, and T the subset containing the remaining elements of T. The links belonging to T are considered as candidates for improvement and must have their three-piece improvement cost functions f(Z) included in the formulation. Let (Z1, X1, C1, C1, Y1, W1) be the solution of the LP where all nonlinear functions f(Z) in (6) are substituted by aZ + J3 , for each link of T. We obtain the following BLP problem with relatively few lower variables:

11

MINZ ZagS Fa(Za,Wma) + Za£T1 Fa(Za,Wma) + Za£T2 (aQZa+l3a) + xG(C) + (l-x)G(C) such that [Z, C, C, X and Y} and {Wma, m=l,Ma, aeS} solve:

H(X,C,C,X,Y) > 0 Z, C, C, X, Y > 0

H5a(Za>Wma) * °

Wma * ° (7)

and such that {Wma, m=l,Ma, acT1} solve:

MIN EaeTl Zm=l,Ma Wma subject to:

Wma " Za - "^ma Wma * °-

where Ma is the number of breakpoints in the piecewise improvement cost function of link a; Fa is the upper objective function of nonconvex improvement cost on link a; and, Hca gives the corresponding constraints.

Let us call (Z2, X2, C2, C2, Y2, W2) the solution of (7). If there exists a

1 ") link j such that j does not belong to T and Z- is different from 0, then j

has to be subtracted from T , added to T , and (7) has to be resolved until all

? 1

links a having Za positive are elements of T , in which case the exact

solution of (6) is found.

We have constructed from the main BLP problem (2) two other secondary BLP formulations: a parameterized BLP (6) and a BLP with few lower variables (7). Problem (7) is solved to give a feasible and optimal solution to (6). If this solution is user-optimized (i.e. solves the lower LP problem (5) when Z is fixed at its optimal value for (7)), then it is feasible for (2) and a higher value of t can be tried. Otherwise, we continue our search for the optimum by decreasing the value of t.

For any solution (Z1, X1, C , C , Y1, W1), we have to distinguish between two different values. The first value is that of the upper objective function

12 of formulation (2), equal to F(Z,W) + G(C); let us call it 01. The second value is the actual cost incurred by the system given a specific value Z1 of Z; let us call it A1. A1 includes actual improvement cost and actual system travel cost. The actual improvement cost is obtained by plugging Z1 in the actual improvement functions f(Z), whether they are linear, convex, or nonconvex. The actual travel cost, on the other side, is obtained after we solve the lower LP problem (5) with Z fixed at Z1, obtain the solution (XJ, C J , C J , Y J ) and plug the XJ into the convex system travel costs; in other words, A1 is equal to fCZ1) + G(CJ). The first value 01 corresponds to a potential optimal, but not necessarily feasible solution; it allows to find ideal solutions to the BLP problem (2). The second value A1 corresponds to a feasible, but not necessarily optimal, solution; it allows to find incumbent solutions for (2).

The proposed algorithm is an iterative procedure that tries at each itera- tion to reduce the gap between ideal solution and incumbent solution. The procedure terminates when the gap is brought below a desired value e, or when the number of iterations exceeds a fixed limit N. The following flow chart gives a general description of the iterative procedure; more details are provided by Appendix B:

Searching for x (Using Binary Search)

	Start with t=1		-
	A. Solve (6)
	B. Test whether solution of (6) is feasible to (2)
	Feaj	sible	Inf easib!	Le
	r 1 Mn1.n t 1«-~„_		PO M«1, ~ n- _~,.11M.~
	Li 1 .	i iom	* i iai 6ej	.	\j£. .		>.e	oiua iici

13 The most difficult step in the algorithm is step A, in which (6) has to be solved. Using the concepts explained above in this section, the solution of (6) is obtained as shown in the following flow chart:

t Fixed

	Al. Start with TX=0, T2=T
	A2. Solve (7)
	A3. Test whether new elements should be added to T
New	elements added	No element?	; added
	A/. M~r» Tl r^r.T T2		K ( C\ r.~1,r~^
			J. , new

Several improvements can be added to the algorithm to make it more effi- cient. The starting ideal and incumbent solutions can be found in a more effi- cient way. A better ideal solution for (2) can be obtained when solving the LP obtained after substituting at x = l all nonlinear functions f(Z) in (6) by ccZ+B; there is no situation better for the system than having cheap improvement and system-optimized flow. Also, the actual cost incurred by the system when the problem is solved as a BLP must be lower than the actual cost incurred when a simpler technique is used; otherwise, there is no justification for using complex techniques. Therefore, to obtain an incumbent solution, we can just solve the problem as an LP; we approximate improvement cost functions by their best linear fit, ignore the lower objective and solve the problem. Let us call (Z°, X°, C°, C°, Y°) the obtained solution, A0 is the value of the starting incumbent solution. Using those techniques makes the gap between starting ideal

14 solution and starting incumbent solution as low as 3 % instead of the 111 % found earlier in the previous section.

Moreover, since the procedure is iterative, it is more efficient to use the information from the previous iterations in order to make the next one easier to solve. One way to do this is to introduce to the problem the constraint:

L < F(Z,W) + G(C) < U

where L and U are the values of the current ideal and incumbent solutions, respectively.

Similarly, the solution of step 2 may require solving (7) more than once.

9 9 9 9 9

Let us assume that we solved (7), obtained the solution (Z , X , C , C , Y ,

9

W ) and have to add a new element and resolve (7). We use the fact that the new

cost will be at most as low as the value we got for the objective when we used

9 1

cheap linear improvement costs (a-Z- + J3-) for the links j to be added to T ;

which gives a lower bound for the solution. An upper bound can also be found by

9

noticing that the new cost is at least as high as the actual cost f(Z )

9

incurred by the system for an added capacity equal to Z . The difference

9

between upper and lower bounds, when step 2 has to be repeated, is E^ti (ai^i + Bj) - f(Z2).

Finally, a feasible solution to a BLP problem is obtained by finding the optimal solution to the lower problem after substituting values for the upper variables. In our case, any optimal solution (Z1, X1, C1, C , Y1, W1) of the BLP problem (6) satisfies the optimality requirement of W1 to the lower problem (A); therefore (Z1, X1, C1, C1, Y1, W1) is a feasible solution of the BLP problem (2) if it satisfies the optimality requirement for the second problem (5). Even when the solution (Z1, X1, C1, C1, Y , W1) is not feasible a new solution (Zk, Xk, Ck, Ck, Yk, Wk) can be obtained, with feasibility guaranteed.

15 It is known from the theory of Linear Programming [see for instance Chvatal 1983, Chapter 10 pp. 148-168] that given two LP problems LP1: (MIN ex, Ax = b, x > 0} and LP2: (MIN ex, Ax = b' , x > 0), where b and b' are different, any basic feasible solution x of LP2 that has the same nonbasic variables as the optimal solution x of LP1, is an optimal solution of LP2. Therefore to generate feasible solutions for BLP problem (2),. which is equivalent to generating optimal solutions to the lower problem, we can solve BLP problem (7) with the additional requirement that all nonbasic variables of the solution (XJ, CJ , C J , Y J ) of the lower problem (5) with Z fixed at Z1 and all nonbasic slack and surplus variables be locked. This is added to the algorithm as shown in Appendix B.

Regarding the reliability of the algorithm, there are no heuristics involved in finding the incumbent and the ideal solutions. The values given by the algorithm for those solutions are exact. As long as the algorithm converges, the solution it gives is the actual optimum within the accuracy desired. However, it is possible that one finds bad examples for this algorithm where the convergence is not satisfactory; in such a case the algorithm ends with a large gap between upper and lower bounds, but does not give a wrong solution.

3. Solution of the Real World Problem

When trying to solve the BLP problem (7) we faced a major constraint, since we could find no reliable BLP software that could handle problems this large. Fortunately, (7) can be transformed into an equivalent Mixed Integer Program- ming problem (MIP), with the number of zero-one variables in the latter equal to the number of lower variables in (7). It should be emphasized that the use of MIP does not mean that it is easier to solve. MIP and BLP are both NP-hard problems and are equally difficult to solve; however, MIP is older and benefits

16 from the availability of good software to solve it. One such package is the Zero/One Optimization Method (ZOOM), which is an extension of the XMP Linear Programming system. ZOOM and XMP are both written in Fortran 77 and revised in 1986 by R. E. Marsten from the Department of Management and Information Systems at the University of Arizona. A description of the XMP library is given in Marsten [1981]. The subroutines of ZOOM and XMP are used by our program whenever an LP or an MIP needs to be solved. The computational study, to be presented in this section, was conducted on the CRAY X-MP/24 supercomputer of the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign.

First, we show how to transform BLP (7) into MIP. The constraint:

Wm = MAX{(Z - qm), 0)

is equivalent to:

either: Wm = Z - qm, or: Wm = 0

which can be equivalently restated, after introducing zero-one variables Vma, as :

Wm - Z < -qm + MVm, -Wm + Z < qm + MVm

(8)

Wm * "(l-V, "Wm * ^-V

where M is a large number (it is sufficient to make M larger than the maximum possible road improvement).

The formulation of (7) as MIP is obtained by replacing the lower problem in (7) by the constraints (8) for each link a. The constraints Wm < M(l-Vm) are trivial because of the nonnegativity of Wm; we take them out, and obtain the following problem to be solved whenever we need to solve (7):

17

«INZ ZaeS Fa(Za'Wma) + EaeTl Fa^a'Wma) + ZaeT2 (°aZa+6a) + *G^ + (l-t)G(C) such that (Z, C, C, X and Y} and {Wma, m=l,Ma, aeS) solve:

H(X,C,C,X,Y) > 0

Z, C, C, X, Y > 0

H5a(Za>Wma) * 0

and such that {Wma, m=l,Ma, aeT ) solve:

(9)

"Wma + Za " MVma * ^ma

Wma " Za " MVma * "^ma

*ma + MVma * M Vma e (0, 1).

The technique used to transform (7) into MIP is different from the technique used by Fortuny-Amat and McCarl [1981]. The use of this latter would have doubled the number of zero-one variables (requiring a zero-one variable for each lower variable and a zero-one variable for each lower constraint).

The large size of the problem (2) is mainly due to the high number (2,095) of routes included in the formulation. Since not all those routes are going to be used, the problem can be simplified by keeping the most attractive ones and discarding the rest. We will solve the problem with the reduced number of routes (we refer to this problem as the small version of the problem), and we solve also the original problem with all routes included; next we will compare the solutions of the two problems. The attractive routes are found by solving both the user equilibrium problem (5) and the system equilibrium problem (same as (5) except that G(C) is replaced by G(C)) as we increase the values of the entries of the trip matrix to make the users and the system choose their next desirable routes. When solving the system equilibrium, we vary also the improvement cost because of its effect on the routes choice. The accumulated number of routes chosen by both the system and the user reaches 389 as the flow

18 fluctuates between the original trip matrix and three times that matrix. Those 389 routes, averaging more than three routes per origin destination pair are the ones to be included in the small version of the problem as compared to the 2,095 included in the large version.

About 15.25 minutes CPU time and 1.4 million words (64 bits per word) computing storage space were required to obtain the solution of the small version of the problem. At iteration 0, we compute starting ideal and incumbent solutions using the ordinary techniques of BLP; the obtained gap is 111.36 %; next, we use the techniques proposed in Section 2 to get better starting solutions. At the same iteration, we solve the problem as an LP and we use the actual cost incurred by the system (namely 257,037.230) as a new incumbent solution. Then we start the iteration 1 with T=l, TflAX13-'- anc* TMIN=^' When we find, in the first step (A1-A2 in the flow chart) of iteration 1, the list of links that are candidates for improvement, we use the solution as a new better ideal solution; the gap decreases to 3.03 %. In the second step (A3-A4-A2 in the flow chart), the MIP (9) is solved with only six integer variables. The solution obtained shows that no more integer variables need to be added; a slightly improved ideal solution is found for the problem (new gap equals 2.93 %) . The feasibility of the solution is tested in step 3 (B in the flow chart); the flow is not user-optimized. We generate a feasible solution using the solution obtained when testing feasibility; a new incumbent solution is obtained and the gap is reduced to 2.57 % at the end of the first iteration. At iterations 2 and 3 there is no success in improving the incumbent solution. Iteration 4, where x = .125, i^x = -250 and t^IN = 0, gives another slight improvement at the end of step 3, when generating a new feasible solution. The procedure is stopped at Iteration 5 because the solution seems to be reaching a steady state; Iteration 5 gives exactly the same solutions as Iteration 4 in all steps. The final solution retained incurs a cost to the system equal to

19 256115.006.

The small version of the problem is not hard to solve even on an ordinary computer. However, the large version is more demanding; it required about 4 minutes CPU time and 2.7 million words computing storage space for running one single iteration. No more iterations were needed because the solution, at the end of the first iteration, was the same as the solution obtained by the small version of the problem at that point; which suggests that our choice for the attractive routes was successful.

Table 1 gives a summary of the intermediate and final results related to the solution of the real world problem; the complete presentation of those results is provided by Appendix E in Ben-Ayed [1988]. The Table includes only the iterations and the steps resulting in an improvement of the gap (increase of the ideal solution or decrease of the incumbent solution). It can be seen from the table that, as is usually the case for iterative algorithms, most progress is made during the first iterations while convergence is slower as the number of iterations increases. The first two columns of the Table are very different from the others because of the instability of the "maximum" variables W in the absence of the corresponding lower objective function; the negative value of the improvement in the first column is a consequence of this instability. The fourth row of the Table specifies whether a solution is ideal or incumbent; when a solution is incumbent, the rows 5, 6 and 7 refer to the actual cost incurred by the system (A1), and when the solution is ideal, those rows refer just to the value of the objective function (01). Starting from the third column, the solution has a tendency to stabilize; the numbers of used links, used routes and improved links are almost the same in all six last columns.

Several results can be inferred from the computational experience and the solution of the real world problem. Concerning the tradeoff between efficiency and accuracy, the equivalence of the solutions of the large version and the

20 small version shows that a complex problem can be simplified without affecting the optimality of the solution. However the simplification must be done carefully and efficiently; other simplifications, such as reducing the problem to a single-level could lead to bad results.

The actual cost incurred by the system when the problem is solved as an LP exceeds the actual cost incurred when a BLP solution is used by 922.224 thousands of Tunisian Dinars. There are two reasons for not having this difference as high as it is supposed to be. First, the travel cost at stable flow is a one-piece linear function that does not depend on the added capacity. However, if stable flow, which is the overwhelming proportion of the flow, does not depend on the added capacity, which is the only instrument in the model for the system to influence the user's choice of the routes, the ability of the system to affect the user-optimized equilibrium is almost paralyzed. Second, the limitation of the added capacity resulted in low improvement costs. It can be seen from Table 1 that improvement costs are about fifty times smaller than travel costs. The introduction of the sophisticated BLP formulation of the improvement costs did not give much saving as compared to the simpler LP formulation because those costs are already low according to the data of the problem. The ability of BLP to give better results is conditioned by the availability of realistic data; providing more representative data to our problem would definitely lead to a larger difference between the solution of the problem when solved as a BLP and its solution when solved as a single level .

Another harmful simplification would be the inclusion of strict capacities on the links, even when the overflow is very small on all links, as is the case in our solution. Imposing strict capacities results in the violation of Wardrop's first condition and therefore the collapse of the equilibrium, which is one of the bases of the problem.

Baene

MEDITERRANEAN

Cape Bon

Pwntslkna (ITALY)

boh Peleg/t (fTALYI

MEDITERRANEAN SEA

100 K.loT^irri

Figure 1 The Links to Bo Improved

22

Summary of the Results
Iteration	0	0	0	1	1	1	1	4 & 5
Step	-	-	-	1	2	3	4	4
Value of t	1.000	1.000	1.000	1.000	1.000	1.000	1.000	< .125
Solution	ideal	incumb.	incumb.	ideal	ideal	incumb.	incumb .	incumb.
Travel Cost	255328.	273736.	252428.	245152.	245320.	252577.	251724.	251680.
Improvement	-124325	3153.	4609.	4329.	4410.	4410.	4424.	4435.
Total Cost	131003.	276889.	257037.	249481.	249730.	256987.	256148.	256115.
Routes Used	278	279	276	275	276	276	275	275
Links Used	99	100	97	97	97	97	97	97
Linkslmproved	18	18	14	14	14	14	14	15
Gap (%)	-	111.36	96.21	3.03	2.93	2.91	2.57	2.56

Table 1

In the final solution, 275 routes and 97 links are used and 15 links need to be improved. Among the 112 links included in the empirical study, 15 are not used by inter-regional traffic; one cannot tell whether those links are redundant in the inter-regional network or their lack of use is due to low values of the corresponding entries in the trip matrix. The study is more concerned with the links to be improved; those are shown on the map in Figure 1. Almost all improvements happen to be on the roads connecting Tunis, the capital, to its major neighboring cities. This result can be intuitively predicted by a simple look at the trip matrix [Ben-Ayed, Blair and Boyce 1988]; the trips originating or ending in Tunis are 56 % of the sum of all the trips of the matrix. A similar recommendation was made by the Plan Directeur Routier [ SETEC 1982 pp. 6-8] stating that the roads corresponding to the exists of Tunis will be highly congested by the year 2000 if no improvements are made. In contrast, the entries of the trip matrix corresponding to regions such as the

23

southern ones are very low, which results in no improvement even for unpaved roads .

		The	Links to B<	2 Improved
Link	Flow	Existing	Added	Improvement	Limit on Added
		Capacity	Capacity	Cost	Capacity
1	2144	1370	774	592	887
3	2537	1637	899	394	899
4	2112	1115	997	297	1083
6	2015	1055	960	615	1070
10	1980	814	1166	484	1460
12	1052	448	604	84	2046
14	306	271	35	22	1354
19	1056	447	609	89	2046
20	2275	1173	1101	463	1101
26	2537	1637	899	245	899
28	2010	1370	640	362	887
30	1517	1370	147	60	887
73	2112	862	1250	394	1496
77	1772	1055	717	156	1070
90	2058	1370	687	179	887

Table 2

In considering the specific link improvement, Table 2 gives for each link to be improved the interregional flow, the existing capacity available to this flow (60 % of the total existing capacity), the added capacity, its cost and the maximum added capacity allowed by the formulation. The improvement required for the link connecting nodes 4 and 29 is very small; a slight improvement of the surface of this link gives it exactly the same capacity as the link connecting nodes 29 and 28 which does not need improvement although it is closer to the capital. All other improvements are more significant and three of them, namely those on links 3, 20 and 26, require the maximum added capacity allowed by the formulation. In each of those links the flow is higher than the new capacity; more detailed study is needed to include the possibility of upgrading them to three or four lane highways. Link 90 is also congested; however the decision maker responsible for the system does not add more

24 capacity although he has the possibility to do so. The reason is that the congestion on those roads is less expensive than the addition of more capacity.

4. Conclusions

BLP remains a little known technique, rarely used despite all its potential capabilities. The computational difficulties of BLP may not be the main reason of its infrequent applications; there are other optimization techniques, such as Mixed Integer Programming, that are not easier to solve but are widely used. BLP is a new model; a great effort is needed to help people know it and take advantage of its powerful formulations. Developing efficient and reliable software that uses data structure and that includes several heuristics and branch and bound algorithms would have considerable effect on increasing the use of BLP.

The empirical study conducted in this research demonstrated that BLP is tractable and can be efficiently used to formulate and solve even complex real world problems such as the design and the improvement of a transportation network. The complexity of BLP should not be a reason for formulating bilevel problems as single- levels ; the solution to the single-level problem might happen to be far away from the actual optimum, and the gain in simplicity and computation time may not be worth the loss in optimality. By taking into account the special structure of the problem at hand and making use of the powerful computers available to us nowadays, it is always possible to solve the BLP problem and get a better solution than the one obtained when approximating it by a single level. For any BLP problem, the use of any related algorithm always gives a solution that is better than, or at worst the same as, the solution obtained when ignoring the lower objective and solving the problem as a single-level. In other words, the relative inefficiency of BLP algorithms does not justify the formulation of bilevel problems as single-levels.

25 Appendix

A. Definition of the Variables and the Functions in Formulation (2)

F(Z,W) = Ea=liU2 Fa(Za>Wma) = Ia=ljll2 ((blaZa+b0a) + Em=1>Ma(bm+lja-bma)Wma

F(W) = ^a=ljii2 Ia(Wma) = Za=l,112 Zm=l,Ma Wma

W) = Za=l,112 GaCCa) G(C) = Ea=l,U2 Ga(Ca)

H1(Y) = 0 is defined as: H2(X,Y) = 0 is defined as: H3(Z,C,X) > 0 is defined as:

= Za=l,112 ^Ca--4c0aka) = Za=l,112 ^a

H4(Z,C,X) > 0 is defined as:

H5(Z,W) > 0 is defined as:

H^(Z) > 0 is defined as:

for each origin-destination pair od:

ZreRod Yr ~ uod for each link a:

Xa " Zr=l,2095 5arYr ~ ° for each link a:

Ca " <WXa + -4ka) * °

ca " cla(xa + -4ka) + (cla " <Wza + dla * ° for each link a:

Ca * c0a<xa + -4ka)

Ca * ^la^a + -4ka) " tela " ^Oa)Za + <ha for each link a:

Wma " Za + %,a ^ ° m=1>Ma

for each link a:

Z„ ■- q3r > 0

'a ^3a Z, C, C, W, X, Y > 0 is defined as: for each link a:

Za> Xa' Ca' ^a> Wma - °

m=l,M.

for all reR

od:

Yr > 0

where Cq3 and c^a are the slopes of the system travel cost on link a at stable and unstable flows, respectively; cQa and cla are the slopes of the cumulative user travel cost at stable and unstable flows, respectively; dla and dla are

26 the intercepts of the travel cost at unstable flow of the system and the cumulative user, respectively; Mfl is the number of breakpoints in the improve- ment cost function in link a; Mfl is 0 when the curve is linear, 1 when the curve has 2 pieces and 2 when it has 3 pieces; bga is the intercept of the first piece of the improvement cost function on link a; bma is the slope of the piece delimited by qm_^ a and qma, m=l,Ma+l, qQa being equal to zero and q^a+^ equal to q3a, the maximum improvement allowed by the model; kfl is the existing capacity of link a; Xa is the flow on link a; Cfl is the system total travel cost on link a; C_a is the cumulative user total travel cost on link a; Za is the number of passenger car units (PCU) added to the capacity of link a; W_a is the maximum of (Za~qma) and 0; Y is the flow on route r; R j is the set of all routes from origin o to destination d; uQ(^ are the entries of the trip matrix; 6ar is binary number equal to 1 when link a belongs to route r, and equal to 0 otherwise.

B. Description of the Algorithm

When describing the algorithm, we use T, T and T , 01 as defined in Section 2:

Iteration 0

• Find a current ideal solution, a current incumbent solution and calculate the current gap for BLP problem (2).

• Call Termination test.

• Set: i = 1, xMAX = 1, tmin =0, 1=1.

• Start iteration I. Iteration I

Set: T1 = <f>, T2 = T. Step 1

• Solve the LP obtained after substituting all nonlinear improvement functions by aZ + 6 for all links to find the links of T that are

27 candidates for improvement.

• Get the solution (Z1, X1, C1, C1, Y1, W1).

• For every link j of T , if the corresponding Zi is not equal to zero,

7 1

subtract j from T and add it to T .

Step 2

• Solve the BLP problem (7) with the set of lower variables equal to {Wma:

bm+l,a-bma < °> m=1>25 asTl)'

• Get the solution (Z2, X2, C2, C2, Y2, W2).

• If there are elements j of T such that Z- is positive, then:

9 i

Subtract from T all elements j and add them to Tx Repeat Step 2 Else:

The optimal solution of (6) is obtained

End If.

2

• If x = 1 and 0 > value of current ideal solution, then:

Update ideal solution and calculate new gap Call termination test End If. Step 3

• Solve the lower LP problem (5) with Z equal to Z.

• Get the solution (X3, C3, C3, Y3) .

9

• If A < value of current incumbent solution, then:

Update incumbent solutions and calculate new gap

Call termination test End If.

• If A2 = 02 (or A2 ~ 02), then:

(Z2, X2, C2, C2, Y2, W2) is feasible

Set: x

MIN " T> x = (TMAX + T)/2

28 Else:

Try to improve feasible solution Set: xMAX = T, t = (xMIN + t)/2 End If. •1 = 1 + 1. • Start iteration I. Termination Test

• If I > N, then:

Number of iterations exceeded before reaching desired accuracy

Incumbent is the best solution found

Stop Else if gap < e:

Incumbent solution is optimal within the specified accuracy e

Stop End If.

• Return.

Try to improve feasible solution

• Solve (6) while locking the nonbasic variables (including slack and surplus variables) of the LP problem (5) solved with Z = Z .

• Get the solution (Z4, X4, C4, C4, Y4 , W4).

• If 0 < value of current incumbent solution, then:

Update incumbent solution and calculate new gap Call termination test End If.

References

BARD, J. F. 1983. An Efficient Point Algorithm for a Linear Two-Stage Optimiza- tion Problem. Operations Research 31, 670-684.

29

BARD, J. F. and J. E. FALK 1982. An Explicit Solution to the Multi-Level Programming Problem. Computers and Operations Research, 9, 1, 77-100.

BARD, J. F. and J. T. MOORE 1987. A Branch and Bound Algorithm for the Bilevel Programming Problem. Department of Mechanical Engineering, University of Texas, Austin, TX.

BEN-AYED, 0. 1988. Bilevel Linear Programming: Analysis and Application to the Network Design Problem. Ph.D. Thesis. Department of Business Administra- tion. University of Illinois, Urbana-Champaign.

BEN-AYED 0. and C. E. BLAIR 1988. Computational Difficulties of Bilevel Linear Programming. Submitted to Operations Research.

BEN-AYED 0., D. E. B0YCE and C. E. BLAIR 1988. A General Bilevel Linear Programming Formulation of the Network Design Problem. Transportation Research B22.

BEN-AYED 0., C. E. BLAIR and D. E. B0YCE 1988. Construction of a Real World Bilevel Linear Program of the Highway Network Design Problem. Submitted to Transportation Research.

BIALAS, W. F. and M. H. KARWAN 1982. On Two-Level Optimization. IEEE Transac- tion on Automatic Control AC-27, 1, 211-214.

BIALAS, W. F. and M. H. KARWAN 1984. Two-Level Linear Programming. Management Science 30, 8, 1004-1020.

BIALAS, W. F. , M. H. KARWAN and J. P. SHAW 1980. A Parametric Complementary Pivot Approach for Two-Level Linear Programming. Research Report No. 80-2, Operation Research Program, Department of Industrial Engineering, State University of Buffalo.

BRACKEN J. and J. T. McGILL 1973. Mathematical Programs with Optimization Problems in the Constraints. Operations Research 21, 37-45.

BRACKEN J. and J. T. McGILL 1974a. Defense Applications of Mathematical Programs with Optimization Problems in the Constraints. Operations Research 22, 1086-1096.

BRACKEN J. and J. T. McGILL 1974b. A Method for Solving Mathematical Programs with Nonlinear Programs in the Constraints. Operations Research 22, 1097- 1101.

CANDLER W. and R. NORTON 1977. Multilevel Programming and Development Policy. World Bank Staff Working Paper 258, IBRD, Washington, DC.

CANDLER, W. and R. T0WNSLEY 1982. A Linear Two-Level Programming Problem. Computers and Operations Research 9, 1, 59-76.

CANDLER, W. and R. TOWNSLEY 1985. Optimizing Subsidy Payments in Linear Programming. IEEE Proceedings of the International Conference on Cyber- netics and Society, 930-934.

CASSIDY, R. G., M. J. L. KIRBY and W. M. RAIKE 1971. Efficient Distribution of Resources Through Three-Levels of Government. Management Science 17, 8,

30

B-462-B-473.

CHVATAL, V. 1983. Linear Programming, Freeman, New York.

DESILVA, A. H. 1978. Sensitivity Formulas For Nonlinear Factorable Programming and Their Application to the Solution of an implicitly Defined Optimiza- tion Model of U.S. Crude Oil Production. D.Sc. Thesis, George Washington University, Washington, DC.

FALK, J. E. 1973. A Linear Max-Min Problem. Mathematical Programming 5, 169- 188.

FALK, J. E. 1977. A Nonconvex Max-Min Problem. Naval Research Logistics Quarterly 24, 3, 441-450.

FALK, J. E. and G. P. McCORMICK 1982. Mathematical Structure of the Interna- tional Coal Trade Model. U.S. Department of Energy Report DOE/NBB-0025 , Washington DC.

FORTUNY-AMAT, J. and B. McCARL 1981. A Representation and Economic Interpreta- tion of a Two-Level Programming Problem. Journal of Operational Research Society 32, 783-792.

GALLO, G. and A. ULKUCU 1977. Bilinear Programming: An Exact Algorithm. Mathematical Programming 12, 173-194.

LEBLANC, L. J. and D. E. BOYCE 1986. A Bilevel Programming Algorithm for Exact Solution of the Network Design Problem with User-Optimal Flows. Transpor- tation Research 20B, 3, 259-265.

MARCOTTE, P. 1986. Network Design Problem with Congestion Effects: A Case of Bilevel Programming. Mathematical Programming 34, 142-162.

MARSTEN, R. E. 1981. The Design of the XMP Linear Programming Library. ACM Transactions on Mathematical Software 7, 4, 481-497.

PAPAVASSILOPOULOS, G. P. 1982. Algorithms for Static Stackelberg Games with Linear Costs and Polyhedra Constraints. Proceedings of the 21st IEEE Conference on Decision and Control, 2 and 3, 647-652.

SETEC Economie et SOTINFOR 1982. Plan Directeur Routier. Rapport General Tomes 1 et 2, Republique Tunisienne, Ministere de 1 ' Equipement , Direction des Ponts et Chausees, Sous-Direction des Etudes.

SHERE, K. and J. WINGATE 1976. Allocation of Resources to Offensive Strategic Weapon Systems. Naval Research Logistics Quarterly 23, 1, 85-93.

WARDROP J. G. 1952. Some Theoretical Aspects of Road Traffic Research. Proceed- ings , Institution of Civil Engineering, Part II, 1, 325-378.

ECKMAN

INDERV INC.

JUN95