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
Live Music Archive
Librivox Free Audio
Metropolitan Museum
Cleveland Museum of Art
Internet Arcade
Console Living Room
Books to Borrow
Open Library
TV News
Understanding 9/11