DECOMPOSITION TECHNIQUES FOR
THE TRAFFIC ASSIGNMENT PROBLEM
By
SIRIPHONG LAWPHONGPANICH
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA IN PARTIAL
FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1983
ACKNOWLEDGMENTS
The author wishes to express his gratitude and appreciation to the
members of his supervisory committee, Dr. H. P. Benson, Dr. D. J.
Elzinga, Dr. R. L. Francis, Dr. H. W. Hamacher, Dr. D. W. Hearn, Dr. R.
W. Swain, and Dr. S. Tufekci, for their assistance and guidance.
The author is especially grateful to his advisor, Dr. D. W. Hearn,
for suggesting the research topic and for his encouragement and
patience throughout the research.
This research was partially supported by the National Science
Foundation under NSF Grant ECS7925065.
TABLE OF CONTENTS
Page
ACKNOWLEDGEMENTS . . . . .. .
ABSTRACT . . . . . . .
CHAPTER
1 INTRODUCTION AND OVERVIEW . . . . 1
1.1 Introduction . . . .
1.2 Overview . . . .
2 PROBLEM FORMULATION AND LITERATURE SURVEY .
2.1 Problem Formulation . . .
2.2 Literature Survey . . .
. 1 . 5
. . 11
3 TRANSFER DECOMPOSITION . . . .. .18
3.1 Partitioning of the Feasible Region in
Transfer Decomposition . . .
3.2 Validation of Transfer Decomposition for the
3.3 Properties of the Mapping n(z) . .
3.4 Algorithms for the Master VIP . .
3.5 Applications: Link Abstraction and SubArea
Focusing . . . . .
. 18
VIP 22
. 28
. 38
. 48
4 APPLICATIONS OF SIMPLICIAL DECOMPOSITION . ... 58
4.1 The Restricted Simplicial Algorithm .
4.2 Simplicial Decomposition for the Asymmetric
Problem . . . . .
. 61
. 78
5 SUMMARY AND FUTURE RESEARCH . . . ... 95
5.1 Summary of Results . . . .
5.2 Summary of Applications . . .
5.3 Future Research . . . .
REFERENCES . . . . . .
BIOGRAPHICAL SKETCH . . . . ... .. ... .107
. .102
2
. . 1
. . 2
. 5
Abstract of Dissertation Presented to the Graduate Council
in Partial Fulfillment of the Requirements
for the Degree of Doctor of Philosophy
DECOMPOSITION TECHNIQUES FOR
THE TRAFFIC ASSIGNMENT PROBLEM
By
SIRIPHONG LAWPHONGPANICH
August, 1983
Chairman: Dr. Donald W. Hearn
CoChairman: Dr. D. Jack Elzinga
Major Department: Industrial and Systems Engineering
This study is concerned with decomposition techniques for the
traffic assignment problem. In particular, it generalizes earlier
decomposition results for the standard traffic assignment problem in
which the travel cost function has a symmetric Jacobian to the case
where the Jacobian is asymmetric. The standard traffic assignment
problem is often formulated as an optimization problem, whereas the
asymmetric case can be formulated in several ways. However, the varia
tional inequality formulation seems to be the most useful and is adopt
ed for this study.
Decomposition ideas originally developed for the optimization
formulation are extended to the variational inequality problem. New
convergent algorithms for both forms of the traffic assignment problem
are given and explored computationally. Particular attention is given
to computer storage requirements and the application of the projection
technique to the master problem of simplicial decomposition. Transfer
decomposition is developed for problems which are only partially
asymmetric. This technique decomposes the problem into two variational
inequality problems, one of which is equivalent to the standard traffic
assignment problem. Furthermore, the impact of transfer decomposition
as a tool in network aggregation is also explored.
CHAPTER ONE
INTRODUCTION AND OVERVIEW
1.1 Introduction
The traffic assignment problem (TAP) is a multicommodity flow
problem which arises when users (or commodities) must share a common
resource, the street network. In this case, each user has a travel
plan; i.e., he must leave from one location, his origin, to another
location, his destination, in the network. It is generally assumed
that users always choose the best possible route under the prevailing
traffic conditions. However, each user's choice is dependent upon
other users' choices as well, because the travel time on each street
depends on the total number of cars on that street: more congestion
means longer travel times. The problem is to determine the route that
each user will choose to complete his trip, which, in turn, determines
the amount of traffic on each street. In transportation planning,
these numbers serve as a prediction of how the traffic as a whole
distributes itself over the entire network, and based on this predict
ion decisions to expand the current street network are made. Present
ly, the TAP is used as the core of many traffic planning systems.
Although reference to the TAP is used throughout this study,
there are other network models in which congestion plays an important
role: communication systems, railway systems, and military logistic
systems. Thus the discussion in this study also applies to other
2
nonlinear cost, multicommodity flow problems as well. One interesting
feature of TAP, not shared by other models, is that it permits several
equivalent formulations.
Since in practice even a moderately sized problem is too
large for many computers, we are interested in decomposition techniques
for the TAP. However, it should be emphasized that the decision to
decompose a problem should not depend solely on the size of the prob
lem. Decomposition should also be considered when the problem exhibits
special structure and when it is desirable to use a special purpose
computer code to achieve some goal (e.g., high numerical accuracy) at
least for part of the network. Such special codes are often incapable
of solving large problems.
This study considers the extension of decomposition techni
ques for optimization problems to the variational inequality formula
tion of the TAP. As shall be seen, the variational inequality formu
lation can be converted to an optimization problem whose objective
function is not necessarily convex or differentiable. However, this
establishes a relationship between a variational inequality and an opti
mization problem which provides evidence for the possibility of extend
ing decomposition techniques in nonlinear programming literature to the
variational inequality problem. Our main goal is to explore this pos
sibility and to unify the ideas for decomposing the variational inequa
lity formulation of the TAP.
1.2 Overview
Chapter Two discusses different formulations of the TAP and
surveys the past work. From the survey, it is evident that this disser
tation is one of a few systematic studies of decomposition techniques
for variational inequality problems. Chapters Three and Four present
two complementary decomposition schemes: transfer and simplical decom
position, respectively. Transfer decomposition is a resourcedirected
decomposition which can be viewed as a special case of generalized Ben
ders' decomposition [Geoffrion 1972]. It was first developed for the
convex programming formulation of the "standard" TAP. In Chapter 3,
this technique is extended to the variational inequality problem.
After the validation of this decomposition scheme, we investigate pro
perties of the resulting master problem which later lead to a computa
tional algorithm for its solution. In the final section of this chap
ter, an application of transfer decomposition to a problem in transpor
tation planning is discussed.
Chapter Four presents simplicial decomposition, a price
directed scheme originally developed for convex programming problems.
Also, it is a generalization of DantzigWolfe decomposition and has
been investigated by many authors in a nonlinear programming context.
First, we investigate further algorithmic enhancements for the existing
simplicial decomposition scheme with a special concern for large prob
lems. In this restricted version, users may specify the number of grid
points to be retained from one iteration to the next. Based on our
computational results, this allows a tradeoff between computer storage
and the progress of the algorithm. In addition, the modified version
also allows for approximation of the master problem's solution which in
effect provides theoretical justification for the ongoing implementa
tions of simplicial decomposition. In the subsequent section, the
standard simplical method is extended to the variational inequality
problem. It is well known that the algorithm may cycle when column
4
dropping is permitted in every iteration. However, in view of the size
of problems normally considered by traffic planners, column dropping is
essential to any column generating scheme; therefore, we propose a col
umn dropping rule. Under this rule, columns can be dropped whenever
the algorithm makes satisfactory progress toward the solution. As in
the restricted case, only an approximate solution to the master varia
tional inequality is needed to establish convergence.
Finally, areas for possible future research are suggested in
Chapter Five.
CHAPTER TWO
PROBLEM FORMULATION AND LITERATURE SURVEY
2.1 Problem Formulation
The traffic network can be represented by a directed graph,
G(N, L), where N is the set of nodes with cardinality n and L is the
set of links or arcs with cardinality Nodes and arcs usually repre
sent intersections and streets, respectively, and a subset of these
nodes serves as origins and destinations for travel demands. Associ
ated with each arc in the network there is a travel cost function which
depends on the aggregate flow on the network. The object is to assign
the travel demands along paths connecting origins and destinations so
as to satisfy a certain criterion. Two widely used criteria first pro
posed by Wardrop [1952] are the system and user optimized criteria.
The system optimized criterion assumes that there is a central authori
ty which controls the traffic, e.g., in a railroad network; thus, it is
only logical that this authority distributes travel demands so that the
total travel cost or delay for the entire network is minimized. On the
other hand, the user optimized criterion assumes that there is no cen
tral authority, and the users themselves decide the distribution of
flows. Naturally, users always choose the cheapest path to reach their
respective destinations. However, this does not mean that one user's
decision is totally independent of the others; in fact, one user's de
cision may affect the decision of every other user since the travel
6
cost functions depend on the aggregate flow. Under this assumption, an
assignment is said to be a user optimized assignment if it prevents
users from wanting to change their assigned routes.
Since the two criteria do not generally yield the same solu
tion, the selection of criteria depends on the underlying application
of the resulting model. For the traffic network, the solution obtained
from the model is usually used as a prediction of the traffic volume on
the streets. Decisions to expand existing streets or to construct new
streets often depend on this prediction. Since the transportation
planners have no control over users of the network (drivers), the user
optimized criterion is the most suitable choice.
Mathematically, one can formulate the traffic assignment
problem with user optimized criterion as a nonlinear complementarity
problem (NCP) or a variational inequality problem (VIP); they are equi
valent to each other. Define
k k
(i) (s t ) = the kth pair of origin and destination nodes
(OD pair) with nonzero demand, and there are K such
pairs.
k
(ii) f. = the amount of travel demand for the kth OD pair
1
which traverses arc i; in other words, it is the
flow on arc i due to the kth OD pair. Let fk
k k k
denote the vector (fl f2 fk ).
(iii) F = the aggregate (link) flow vector (or flow
pattern) induced by fk, k = 1, ., K;
that is F = (FI, F2, ... F) and F. = f.
k k k
(iv) P = the set of paths from s to t which has
cardinalit k and P = M = k m
cardinality m Let P = U P and P = M m
(v) h
(vi)
(vii)
d
C(F)
(viii) B (B)
= the path flow vector where the jth element of h
represents the flow on the jth path from the set P.
k k
= the (fixed) travel demand from s to t
= the travel cost vector and C(F) =
(C (F),C2(F),..., C (F)).
= the nodearc (arcpath) incidence matrix,
where
B1
Bij =
13 0
I 1
0
if node i is incident to arc j
otherwise
if arc i belongs to path j
otherwise
Then, the vector F is said to be a feasible flow vector if it belongs
to the following set:
k k k k
(2.1) S = { F : F. = f. Bf = r f. > 0, for all i and k }
1 1
k k
where r the requirement vector for the kth OD pair, has two nonzero
entries with values dk and dk Also, the vector h is said to be a
feasible path flow vector if it satisfies the following equations:
(2.2) E h. = dk k = 1,..., K
jePk
(2.3) h. 2 0 V j P.
Obviously,
(2.4) F = B h,
and the travel cost function can be defined in terms of path flows as
follows:
(2.5) (h) = BT C(F)
The travel cost vector has, until recently, been assumed to
have the form C.(F) = C.(F.), i.e., the travel cost on link i depends
only the flow on link i, F.. Under this assumption, the Jacobian of
1
C(F) is diagonal, hence symmetric, and, as explained in the next
section, the TAP is equivalent to a convex optimization problem.
However, in a realistic situation, the travel cost on each link of the
network also depends upon the flow on other links as well. An example
of such a situation arises in uncontrolled traffic intersections where
the travel cost on a link depends on the traffic on all intersecting
links in an arbitrary manner, i.e., the influence of flow on link i on
the travel cost of link j is not necessarily the same as the influence
of flow on link j on the travel cost of link i. Therefore, the
Jacobian of C(F) should be asymmetric in order to have a more realistic
model. This model has been studied both as a nonlinear complementarity
and a variational inequality problem, and they are the subjects of the
next two subsections.
2.1.1 Nonlinear Complementarity Formulation
It can be argued (e.g., Aashtiani and Magnanti [1979]) that a
path flow vector satisfies the user optimized criterion if and only if
all utilized paths between a given OD pair have equal travel cost which
is lower than or equal to the lowest travel cost among the nonutilized
paths. Under this situation, no user has any incentive to take a dif
ferent path, and the user optimized path flow can be viewed as the
steady state or equilibrium path flow vector. Formally, h is a user
optimized or an equilibrium path flow vector if it satisfies
k
(2.6) (C.(h) u ) = 0 if h. > 0 for j c P, k = 1, ..., K,
Sk 3
(2.7) (C.(h ) uk) 0 if h = 0 for j e k = 1, ..., K,
k3
(2.8) E h = d k = 1 ,..., K,
(2.9) h. > 0 V j C P,
k k
where uk is the equilibrium travel cost between s and t Let
x = (h,u) ,
(2.10) y.(x) = C(h) u j P and k = 1, ..., K,
k
(2.11) Zk(x) = Z h. d k = 1, ..., K,
jkpk 3
and
(2.12) G(x) = (Yl(x),...,yM(x),Zl(x),...,ZK(x)).
Then, the equations (2.6) to (2.9) are identical to
(2.13) G(x)x = 0
(2.14) G(x) > 0
(2.15) x 0 ,
which is the standard form of a NCP. (For convenience, we denote the
inner product of two vectors, a and b, by ab, i.e., ab = E a.b..)
i 1
2.1.2 Variational Inequality Formulation
Smith [1979] gave a dynamic interpretation to the user
optimized criterion as follows: Consider two time epochs today and
tomorrow. Users are not permitted to change their routes during today,
but may do so tomorrow provided that the total travel cost, based on
today's path time (cost), is reduced. Note that this interpretation
also permits each user to change his route if he individually would
10
reduce his travel cost, since any such change must reduce the total
*
travel cost. Thus, h is an equilibrium path flow vector if no user
can find a faster route based on today's travel time. Mathematically,
* '.
h must be an element of S = { h : h satisfies equations (2.2) and
(2.3) } and
(2.15) C(h )(h h ) V 0 V h c S
or, in terms of link flows, F must be an element of S and
(2.16) C(F )(F F ) 2 0 V F e S.
*
The problem of finding h (F ) satisfying the conditions stated above
is called a variational inequality problem and is denoted as
VIP [S,C] (VIP[S,C]).
The following theorem from Smith [1979] establishes the
equivalence of the NCP and VIP formulations:
"i".,
Theorem 2.1 The vector h is a solution to VIP[S,C] if and only if
it satisfies equations (2.6) (2.9).
2.1.3 Equivalent Optimization Problems
It is wellknown [Ortega and Rheinboldt 1970] that if C(F)
has a symmetric Jacobian, then C(F) is a gradient mapping for the
following functional:
(2.17) Z(F) = f1 (F F ) C(FO + t(F F )) dt
0
where FO is an arbitrary vector in the domain of C(F), and the NCP and
VIP formulations can be viewed as the KuhnTucker conditions for the
following optimization problem [Dafermos 1971]:
min { Z(F) : F S } .
P2.1:
11
Since (h) = C'(F) B, where the matrices (h) and C'(F) denote the
Jacobian of C(h) and C(F), respectively, C(h) has a symmetric Jacobian
if and only if C(F) has one also. Therefore, an optimization problem
analogous to P2.1 exists for the TAP formulated in pathflow space as
well.
When the Jacobian of C(F) is asymmetric, Hearn and Nguyen
*
[1982] showed that F is an equilibrium flow if and only if F solves
P2.2: min { G(F) : F S } ,
where G(F) = max { C(F)(F H) : H e S }. This function G(F) is
called the gap function in Hearn [1982]; generally, it is neither
convex nor differentiable. By further assuming that C(F) is monotone,
i.e.,
(2.18) (C(F1) C(F2))(F1 F2) Z 0 V Fl, F2 C S,
they also showed that the above result also applies to the following
optimization problem:
P2.3: max { g(F) : F S } ,
where g(F) = min { C(H)(H F) : H E S }. Hearn et al. [1983] show
that G(F) is a dual function of g(F) and vice versa.
2.2 Literature Survey
The survey of decomposition techniques for the TAP is divided
into two subsections, one for the case of symmetric Jacobian and the
other for the asymmetric case. In all of the papers concerning decompo
sition techniques for the symmetric case, the authors assumed that the
Jacobian of C(F) is diagonal, which implies that Z(F) is separable by
arcs. However, the results, in most cases, still hold when the
12
assumption is relaxed. For the asymmetric case, only the VIP formula
tion seems to lend itself to decomposition, since it can be easily
transformed into optimization problems for which there exist many de
composition techniques. Thus, the emphasis here will be on decomposi
tion of the VIP formulation for TAP, a rather new area of research as
far as the transportation literature is concerned. It was only in 1979
that Smith first formulated the TAP as a VIP, and, as a result, there
are only three papers on the subject. In view of the rather small num
ber of papers on this topic, the survey below will include papers on
algorithms for VIP in the transportation literature as well.
In this study, we think of decomposition as a technique which
divides the entire problem into two or more simpler or smaller prob
lems, one of which is often designated as the master problem and the
rest as subproblems. Some methods which are referred to as decomposi
tion in the TAP literature, such as Dafermos [1980a], are better clas
sified as approximation methods, which are discussed below.
2.2.1 The Symmetric Case
Decomposition techniques for the symmetric TAP can be dichot
omized into classes of resource and price directed techniques. The
latter class includes column generation and simplicial decomposition.
Leventhal et al. [1973] applied the standard column generating techni
que to the path flow formulation and obtained a subproblem which gen
erates paths violating the KuhnTucker conditions. Later, Cantor and
Gerla [1974] applied simplicial decomposition to the arc flow formula
tion, i.e., problem P2.1, and obtained as a subproblem a linear cost,
uncapacitated TAP. Furthermore, simplicial decomposition can be viewed
as a generalization of the FrankWolfe algorithm or the DantzigWolfe
13
decomposition, and it is also applicable to general nonlinear program
ming problems with linear constraints. Table 21 summarizes the past
work on this technique.
All of the resource directed techniques are of the general
ized Benders' type, and they partition the underlying network, which,
in this case, is the "resource" to be shared by users. These techni
ques differ from one another in the manner by which the network is
divided. Barton and Hearn [1979] presented two methods for dividing
the network: one is to split it into two or more subnetworks at the
nodes, and the other is to partition the arc set into two or more sub
sets which then define the subnetworks. They called these methods
subset and transfer decomposition, respectively. A special case of
subset decomposition is the geographical decomposition technique [Maier
and Robinson 1977] which first decomposes the problem into K TAPs, each
having only one OD pair with positive demand. Then subset decomposi
tion is applied to each TAP problem with one OD pair, or to use a more
familiar name, a nonlinear minimum cost flow problem. To summarize,
all decomposition techniques in the resource directed class alter the
underlying network structure by decomposing it into two or more sub
networks, and, for this reason, a more appropriate name for them would
be "structural decomposition."
2.2.2 The Asymmetric Case
Bertsekas and Gafni [1982] proposed a projection algorithm
for the path flow formulation of the TAP, and, noting that the number
of paths are too numerous to generate a priori, they devised a scheme
to generate them sequentially. The resulting algorithm is essentially
the simplicial scheme in the path flow space. Concurrent with our
Table 21: Summary of Simplicial
Obj. Function
Feasible Region
Master Problem
Monotonicity
Argument
Column Dropping
Argument
# of major iter.
von Hohenbalken Holloway Sacher Cantor & Gerla Leventhal, et al.
[1975] & [1977] [1976) [1980] [19741 [1973)
pseudoconvex convex quadratic Klienrock's convex
(convex) delay function
polytope compact and polyhedron network network
convex set constraints constraints
nonlinear prog. nonlinear prog. Linear nonlinear prog. nonlinear prog.
Complementarity
van Moeseke's KT cond. KT cond. KT cond. KT cond.
lemma
follows from follows from Caratheodory Caratheodory Uniqueness of
monotonicity monotonicity Theorem Theorem Lagrange mult.
finite infinite finite infinite finite
Decomposition
15
study, Smith [in press a] independently proposed the simplicial
decomposition scheme in the link flow space and his own algorithm
[Smith in press b] to solve the master VIP. In these papers all gene
rated extreme points are retained until the algorithm terminates,where
as our algorithm emphasizes the elimination of superfluous extreme
points whenever possible. Subsequent to our study, Pang and Yu [1982]
viewed the convex hull of retained extreme points during each iteration
as a sequence of subsets of the feasible region and gave conditions for
global convergence of their method. Currently, the sequence of subsets
generated from our algorithm is the only one known to satisfy their
conditions.
In the transportation literature, there are two major classes
of algorithms: one involves solving the VIP using an approximation of
the mapping C(F) and the other uses existing optimization algorithms to
solve the equivalent optimization problems formulated in Section 2.1.
The prototype algorithm for the first class of algorithm is as follows:
Given Ft e S, let Ft+l solve the following approximate
VIP [S, Ct]: find Ft+l S such that
(2.19) C(Ft+ )(F Ft+) 2 0 V F S ,
where Ct() is an approximation of C(.) at the point F.
Below is a list of possible approximations for the vector function,
C(F). Among them, the projection method, originally introduced by
Dafermos [1980b], is the only globally convergent algorithm; the others
only converge locally. All of these approximations are discussed in a
survey paper by Pang and Chan [1982] and a book by Ortega and Rhein
boldt [1970]. Therefore, we only list below the references pertaining
to the TAP.
16
(i) Diagonalization [Abdulaal and Leblanc 1979; Dafermos
1980a, 1982, 1983; and Florian and Spiess 1982]
The function C (*) is defined to be a diagonal function
t t t t t
C. (F) = C. (F,..., F Fi, Fi+1.. F ) i = 1 ..., m
1 1 l' m
and the approximate VIP becomes the standard TAP.
(ii) Linear Approximation [Aashtiani and Magnanti 1979]
The function Ct(*) has the form
Ct(F) = C(Ft) + M(Ft) (F Ft) ,
where M(Ft) is an m x m matrix. Aashtiani and Magnanti let M(Ft) be
the Jacobian of C(F) at F = F
Outside the transportation literature, the choices of M(Ft)
are more numerous, e.g.,
L(Ft) or
M(Ft) = D(Ft) or
(U(Ft) + D(Ft)) / w*
where D(Ft), L(Ft), and U(Ft) are respectively the diagonal, lower and
t *
upper triangular submatrices of the Jacobian of C(F ), and 0 < w < 2.
(iii) Projection Method [Dafermos 1980b, 1983; and Bertsekas and
Gafni 1982]
The projection method is a special case of the linear
approximation method where the matrix M(F ) is taken to be G, some
fixed symmetric and positive definite matrix. The reason for the term
projection is that the vector F solving the approximate VIP [S,Ct] is
precisely the projection of the point Ft G (Ft) onto the set S,
where the projection is defined with respect to the Gnorm, i.e.,
Ft+l = Pr (t G C(F))
S
17
where for a given vector Z, PrS (Z) is the (unique) vector solving the
problem
min { Z Y G : Y S }
where lix I= (xGx).
Chapter 8 of Auslender [1976] and a recent paper by Dupuis
and Nguyen [1982] discuss optimization algorithms to solve the dual
pair of problems presented in Section 2.1. The algorithms discussed
therein include steepest descent, penalty and barrier methods, cutting
plane and subgradient optimization. Also, Nguyen and Dupuis [1981]
improved the cutting plane algorithm originally developed by Zuhovicki
et al. [1969], and applied it to the TAP. Their paper contains the
most complete computational experiments, so to judge the efficiency of
some of our algorithms, we will use their results as a basis for the
comparison.
CHAPTER THREE
TRANSFER DECOMPOSITION
Barton and Hearn [1979] developed the transfer decomposition
technique for the symmetric traffic assignment problem. Also, they
showed that it is a special case of generalized Benders' decomposition
[Geoffrion 1972]. In this chapter, we develop, extend, and generalize
the technique to VIPs over a polyhedral region. When transfer decompo
sition is applied to a general convex programming problem with linear
contraints, the master problem is generally a nondifferentiable convex
programming problem even when the original problem is differentiable.
Analogously, the master VIP resulting from the extension is a general
ized VIP, a harder class of VIPs. However, as in the symmetric case,
it can be shown that the master VIP resulting from the TAP reduces to
the standard VIP for which there exist efficient solution techniques.
This chapter is divided into four sections. The first sec
tion discusses how transfer decomposition partitions the feasible re
gion which is assumed to be in a certain form. The second section
validates transfer decomposition for VIPs. The third section discusses
properties concerning the master VIP. The last section discusses appli
cations of transfer decomposition to two problems of traffic assignment
methodology, subarea focusing and network aggregation.
3.1 Partitioning of the Feasible Region in Transfer Decomposition
Consider the following VIP, VIP[S,(HI,..., H H ,., H )]:
1 p 1 q
^1 ^2 p, ^ 1 ^2 ^q
find (u u ,...,u v v ..., v) S such that
18
19
(3.1)
P q. q
P i i ^i 4 1 i i
E H. (u ) (u u ) + E H. (v) (v v ) a 0
i=l i=l
i=1 1=
V (U1,..., v1,..., vq) E S,
(1 1 V) which satisfies
where S is a set of (u ,...,u ,v v) which satisfies
i i
(i) A.u = E.z i = ,..., p
1 1
P i q i
(ii) Z E.z + E A.v = r
i=l i=l1
(iii) a.. 5 u < b.. ;
11 i
a. v.
13 3 13
where A., A. and E. are
1 1 1
a.., a.., b.., b.., b..
The above VIP
Barton and Hearn. They
problem:
' i '
a.. < z. < b..
13 3 13
V i, j
V i, j ,
matrices, r is a vector of scalars, and a..,
13
are scalars.
includes the class of problems considered by
consider the following multicommodity flow
P3.1 min E Z.(X.) + Z (Y.)
e jcL 2
jEL1 jL2
subject to
k k k
Bx + B2y = r k = 1, ..., K
X k k
k k
k k
x y k 0 k = 1, ..., K ,
where B = [B1, B2] is the nodearc incidence matrix and B1 and B2
contain columns corresponding to arcs in arc sets L1 and L2 which
k k
partition L. The vectors x and yk denote flows on the arcs in L and
L2, respectively. Problem P3.1 can also be stated as a VIP: find a
'Y,
20
feasible point (X Y ) such that
(3.2) E Z.(X *)(X. X.) + Z.(Y )(Y. Y.) 0 V (X,Y) feasible,
jeL13 3 3 jeL2 3 3 3
1 2
where Z.(*) (Z.(*)) denotes the derivative of Z. (*) (Z.(*)). In
terms of the earlier notation, f = (xk y ), F = (X, Y), C(F)
U I
(C (X), C2(Y)), where C1(X) = (Z (X1), Z(X2),...) and C(Y) =
(Z1(Y1), Z2(Y2),...), and S is the set of (X,Y) feasible to problem
P3.1.
To decompose problem P3.1 Barton and Hearn introduce artifi
k
cial variables w k = 1,..., K, and a nodearc incidence matrix D and
rewrite the feasible region as
k k k k k k k
(3.3) S = ( (X,Y) : Dw + B2y = r ; B1x = Dw ; X = x Y = y,
W = E w and x w y 0
Intuitively, the equations in (3.3) describe the process of removing
arcs in L from the network G(N,L) and rejoining the severed network,
G(N N, L \ L ), by artificial arcs designated by w's. Here, N1
denotes the nodes adjacent to arcs in L \ L1. This process transforms
the original network into two networks, one (the subnetwork) consists
of arcs in.L1 and the other (the master network) consists of the remain
ing arcs L \ L and the artificial arcs. Figure 31 depicts an example
of this transfer decomposition process.
In this example,
L1 = { (5,6),(5,7),(5,9),(6,5),(6,8),(6,9), (7,8),(8,7),(9,7),(9,8)}
L \ L1 = { (1,5),(1,6),(2,5),(2,6),(7,3),(7,4),(8,3),(8,4)}
N1 = {1, 2, 3, 4, 5, 6, 7, 8 }
21
Y X Y
(a) Original Network
D

5
3~^t

L2
Master Network
Origins
Destinations
/
x
Subnetwork
(b) Decomposed System
Figure 31: Transfer Decomposition
L 2
L1
L2
L2
22
N2 = {5, 6, 7, 8, 9} = nodes in the subnetwork,
and the set of artificial arcs = {(5,7),(5,8),(6,7),(6,8)}.
Let A1, E and 1 be block diagonal matrices each having B1,
1 1 1
D, and B2, respectively, as its block, the vectors u z and v be
1 K 1 K 1 K
the vector (x ,..., x ), (y ,...,y ) and (w ,...,w ), and finally,
H (u1) and (v1) be the vectors (CI( E xk),..., C( x k)) and
(C2( E yk),...,C2( E yk)). Since
H (u) (u u) + (v v)
^k k ^k ^k k ^k
SE C( x) (x x ) + E C ( y)(y y
k k
= C1(X) (X X) + C2 (Y)(Y Y),
k k
where, as before, X = E x and Y = E y It is easy to establish the
equivalence between VIPs defined by equations (3.2) and (3.1) with
p = q = 1.
3.2 Validation of Transfer Decomposition for the VIP
It suffices to prove the validity of transfer decomposition
as applied to the VIP defined by equation (3.1) with p = q = 1, all
lower bounds equal zero, and all upper bounds equal m. In this case,
the VIP reduces to VIP[S, (H, H)] : find (u ,v ) e S such that
(3.4) H(u*) (u u*) + 1(v) (v v ) b 0 V (u, v) e S,
where
S = { (u,v) : Ez + Av = r, Au = Ez; 0 5 v., u., z. i m}.
Also, it is assumed that
(i) (H(u), f(v)) is continuous, and
(ii) the set
S(z) = { u : Au = Ez O 5 u. 5 m } ? p V (v,z) E S',
where S' = { (v,z) : Ez + Av = r, 0 < vi, z. < m } From (i), there
11
23
must exist at least one solution to VIP[S,(H,i)] (Theorem 4.1, Kinder
lehrer and Stampacchia [1980]).
We decompose VIP [S,(H,i)] as follows:
The master VIP, MVIP [S',(nr)]:
T *
Find (z ,v ) e S' such that for some E T e C(z )
T V *
(E ) (z z) + H(v )(v v ) > 0 V (z,v) C S',
t T *
where n(z ) = {E : H(u ) A ~ X + a = 0, Xu = 0,
*
a.(u. m) = 0 V i, a X. a 0, and u c S (z ) } ,
11 i 1
and u solves the following VIP:
The SubVIP, SVIP[S(z ),H]:
*
Find u e S(z ) such that
*
H(u )(u u ) 2 0 V u e S(z ) .
If H(u) has a symmetric Jacobian, then n(z) is the subdiffer
ential of the function
S*
V(z) = min { h(u) : u e S(z ) } ,
where
h(u) = 1 (u u ) H(u + t(u u ))dt,
'i *
and u would solve the problem V(z ). That is, SVIP[S(z ),H] becomes
an optimization problem, and the above process provides a method which
extracts the optimization component from the original VIP. This could
lead to a significant computational improvement since optimization
problems can often be solved more efficiently than VIPs. The following
theorem (Theorem 2.1, Auslender [1976]) which relates Lagrange multi
pliers (dual variables) to VIPs is essential to the validation of the
above decomposition process.
24
Theorem 3.1 Let S = { u: Au = r, u a 0 } and H(u) be continuous.
Then, u e S satisfies
*
H(u )(u u ) 0 V u e S,
if and only if there exists 5 and A such that
H(u ) A T X = 0
Au = 0
X 0.
Proof: Apply the KuhnTucker Theorem to the following linear pro
gramming problem:
min { H(u )(u u ) : u S } O
The following theorem verifies the above decomposition
process:
Theorem 3.2 VIP[S,(H,H)] and the pair, MVIP[S',(n,H)] and
SVIP[S(z),H], have equivalent solutions.
*
Proof: First, assume that (u ,v ) solves VIP[S,(H,H)]. Then by
*
Theorem 3.1, the pair (u ,v ) must also solve the following linear
programming problem:
m *
P3.2: min { H(u ) (u u ) + H(v )(v v ) : (u,v) S } .
Applying Benders' decomposition to problem P3.2, we obtain the follow
ing pair of problems:
P3.3 min { V(z) + H(v )(v v ) : (z,v) e S' }
and
P3.4 V(z) = min { H(u )(u u) : u e S(z) }
25
*
By Theorem 2.1 of Geoffrion [1972], there exists a z such that (z ,v )
*
solves the master problem, P3.3, and u solves the subproblem, P3.4.
Furthermore, the function V(z) is convex and by Theorem 3 in Section
9.3 of Lasdon [1970] and Theorem 23.9 of Rockafellar [1970], the
subdifferential, 3V(z), is
8V(z) = { E T : H(u) A T A + a = 0; Au = 0; a (u. m) = 0
V i, and a., X. 0 and u solves P3.4 }.
1 1
Thus, both problems, P3.3 and P3.4, are convex programs and their
solutions must satisfy the following optimality conditions: for prob
T *
lem P3.3, there exists a vector E T E v(z ) such that
(t J *
(3.4) (E) )(z z ) + H(v )(v v ) 0 V (z,v) E S',
and for problems P3.4,
(3.5) H(u )(u u ) Z 0 V u S(z ).
*
Since 8V(z ) = n(z ), the equations (3.4) and (3.5) define the master
*
and subvariational inequality problems, respectively. Hence, (z ,v )
and u solves MVIP[S',(n,H)] and SVIP[S(z ),H], respectively.
I',
To show the converse, let (z ,v ) solve MVIP[S',(n,H)], i.e.,
T *
there exists a vector E T c n(z ) such that
T *
(E ) (z z ) + H(v )(v v ) a 0 V (z,v) E S'
or,
,
(3.6) E(Ez Ez ) + H(v )(v v ) 2 0 V (z,v) E S'.
Since the set S(z) < for all (z,v) e S', there must exist at least
*
one u feasible to the set S(z), i.e., Au = Ez and 0 < u. < m. Let u
1
26
** *
solve SVIP[S(z ),H]. Then, substituting Au and Au for Ez and Ez,
respectively, in equation (3.6) gives the following:
S(Au Au ) + H(v )(v v ) L 0
(AT + ( *
(A ()(u u ) + H(v ) (v v ) 2 0
V (u,v) E S
V (u,v) E S.
T *
However, from the definition of n(z ) we have A T = H(u ) X + a
which can be substituted into equation (3.7) to obtain
*
[H(u ) A + a ](u u ) + H(v )(v v ) 2 0 V (u,v) E S
or
*
(3.8) H(u )(u u ) + H(v )(v v ) a (X a )(u u ) V (u,v) E S.
*
Furthermore, (X a )(u u ) = E (X. a.)(u. u.) and
2 1 1
*
(i) if u.
1
*
(X.
(ii) if u.
12
(X.
2.
*
= 0, then a. = 0, X. 2 0 and

 a.)(u. u.) = Lu. 0.
2. 2. 2.2
*
= m, then a. >
*
 a.) (u. ui)
2. 2. 2.
*
0, X. = 0 and
*(u )
= a.(u. m) 0.
2.2.
*
(iii) if 0 < u. < m, then a. = 0, X. = 0 and
1 1 1
*
( a. a.)(u. u.) = 0.
Thus, the right hand side of equation (3.11) is nonnegative and
(u ,v ) solves VIP[S,(H,H)]. O
or
(3.7)
27
3.2.1 An Example
To illustrate the result, consider the following VIP:
S = { (Ul, u2, V, V) : V + = 1, u + u2 = z, Vi, ui, z a0 }
H(u) = (5u1 + 2u2; 3u1 + 6u2)T
rT
H(v) = (3v1+ v2; 2v1 + v2)
The subvariational inequality problem becomes the following:
find u e S(z) = { u: uI + u2 = z; u. i 0 } such that
1
(5u1 + 2u2)(ul ul) + (3u1 + 6u2) (u2 u2) > 0 V u e S(z),
*
and the solution is u = 2z/3 and u2 = z/3. For any z e S', the
mapping n(z) is
n(z) = {A : 4u A + a = 0, 4u X + a2 = 0, au = 0, a. 2 0} = {4u},
and the master VIP can be written explicitly as follows:
find (z,v) e S' = { (z,v) : z + v + v2 = 1; z, vl, v2 0 } such that
4z (z z ) + (3v + v)(v 1) + (2v +.4v2)(v v2) a 0
V (z,v) E S
and the solution is z = 10/26, v1 = 12/26 and v2 = 4/26. This yields
u = 20/78 and u2 = 10/78. Note that
H(u*) = (40/26, 40/26)T and (v) = (40/20, 40/20)T;
H(u ) = (40/26, 40/26) and H(v ) = (40/20, 40/20)
thus, the obtained solution solves the VIP.
28
3.2.2 A Special Case H(u) and H(v) have Symmetric Jacobians
If both H(u) and H(v) are gradient mappings of convex
functions, h(u) and h(v), then n(z) is the subdifferential of the
following convex function:
P3.5: V(z) = min h(u)
subject to
Au = Ez
0 5 u. 5 m,
1
and the master VIP reduces to
P3.6: min V(z) + h(v)
subject to
Ez + Av = r
0 S vi., z. m.
This is exactly what Barton and Hearn [1979] obtain when generalized
Benders' decomposition is applied to the problem below:
P3.7: min h(u) + h(v)
subject to
Ez + Av = r
Au = Ez
0 < u., v., z. < m.
3.3 Properties of the Mapping n(z)
The purpose of this section is to investigate properties of
the mapping n(z). As defined previously, n(z) is generally a pointto
set mapping which further implies that the master VIP belongs to the
29
class of generalized VIPs [Fang and Peterson 1982]. The solution of a
generalized VIP is significantly more difficult than a standard VIP.
We know of only one proposed, but untested, method [Fang and Peterson
1982]. Therefore, it is advantageous to consider conditions under
which the master VIP reduces to an ordinary VIP, i.e., those for which
n(z) is pointtopoint. An important result is that the transfer
decomposition of a TAP always results in a standard VIP for the master
problem.
In addition, we also investigate various monotonicity and
continuity properties of the mapping n(z). These properties are impor
tant in the development of algorithms for the master VIP. The majority
of existing algorithms for VIPs require the mapping to be, at least,
monotone and continuous.
For the sake of completeness, some well known definitions
(e.g., Auslender [1976]) essential to the analysis are listed below.
Definition Let n( *) be a pointtoset mapping defined on S ( R then
(i) n(*) is monotone on S if
(S )(x y) 2 0 V x,y E S and 5 E (x),
x y x
C n(y).
y
(ii) n(') is strictly monotone on S if
(C E )(x y) > 0 V x,y E S and 5 e n(x),
x y x
e n(y).
(iii) n(*) is strongly monotone on S if there is a positive
scalar a such that
30
(5 ) (x y) 2 a(x y) (x y)
x y
V x, y e S and S e n (x), e (y).
If n (*) is a pointtopoint mapping, E and E are replaced by n (x) and
x y
n (y), respectively, in the above definitions.
Theorem 3.3 If H(u) is monotone on the set of all feasible u, then
n(z) is also monotone on the set of all feasible z.
Proof: For all z, = z2 feasible and E TE E n(z ) and E TE e n(z2)
we have
(E TE E T2)(I z2)
= (1 2) (Ez1 Ez2)
= (5 E2) (Au Au2)
T T 1 2
= (AE A 2)(u u)
=(H(u ) 1 + a H(u2) + 12 a2) (u u2)
=(H(u ) H(u) (u u2)
1 1 2 2 2 1
+ (X a 2 a2) (u u1),
1 2
where u and u are solutions to the subvariational inequality prob
1 1 2 2
lems, SVIP[S(z ),H] and SVIP[S(z ),H], and (X a ) and (X a ) are
the Lagrange multipliers associated with the solution vectors u and
2
u respectively. Since H(u) is monotone, the first term of the last
equation is nonnegative. It is shown below that the last term of the
last equation is also nonnegative, and thus n(z) is also monotone.
31
The last term can be written as
1 1 2 2 2 1
Z (X. a A. + a.)(u. u.),
i i i1 1 1 1
1
and we must show that each element of the summation is nonnegative.
There are nine cases to be considered:
1 2 1 1 2 2
(i) 0 < u. < m and 0 < u. < m: Here ai, a., Xi and a.
are zero, thus (1 a .2 + a2)(u u) = 0.
S i 1 1 1) 1
1
(ii) u. = 0 and 0
1
2 2 1
A., and a. are zero, thus (X.
1 1 1
2
< u.
 a1
Here 1
1
2 2
+ a.)(u.
1 1
1
is nonnegative and a.,
1 1 2
 u.) = X.u. > 0.
i i i
1 2 1 1
(iii) u. = m and 0 < u. < m: Here a. is nonnegative and .i,
1 1 1 i
2, and a2 are zero, thus (X1 a X + a2)(u2 u) = a(u m) 0
. i i i i i i 1 11 1
2
(since u. < m).
Similarly, it can be shown that the term
1 1 2 2 2 1
(X. a. A2 + a.)(u. u.) is nonnegative for the following cases:
1 2
(iv) 0 < u. < m and u. = 0.
1 1
(v) u.
(vi) u.
(ii) u
1
1
2
= m and u. = 0.
1
2
= 0 and u. = 0.
1
2
= 0 and u. = m.
1
32
1 2
(viii) u. = m and u. = m.
1 1
1 2
(ix) 0 < u. < m and u. = m. 0
1 1
Corollary 3.4 Assume that S(z ) n S(z2) = ) n(z) is strictly mono
tone if H(u) is strictly monotone.
Proof: From Theorem 3.3,
(ET1 E 2) (z1 z2) 2 (H(u ) H(u 2))(u u2
1 1 2 2 2 1
+ (1 a x + a )(u u )
a (H(u ) H(u 2)) (u u2) > 0.
The last strict inequality follows from the strict monotonicity of H(u)
1 2 1 2
and the fact that u 1 u (since u e S(z ) and u S(z2)).
Corollary 3.5 Assume that E has full column rank. The mapping n(z) is
strictly monotone if H(u) is strictly monotone.
Proof: We show that the sets S(z ) and S(z ) are disjoint for any
z1 f z2. Assume the contrary, i.e., there exist zI and z2 such that
1 2 1 2
z1 i z2 and u = u for some u e S(zl) and u e S(z2). Since
1 2
Au Ezl, and Au = Ez2,
1 2
0 = A(u u ) = E(z1 z2).
By the assumption that E has full column rank, z1 z2 = 0 or zI = z2.
This contradicts the assumption that z1 # z2. Therefore, the sets
S(z ) and S(z2) must be disjoint and the result follows directly from
the previous corollary. 0
Corollary 3.6 Assume that E and A are nonsingular. n(z) is strongly
monotone if H(u) is strongly monotone.
33
Proof: From Theorem 3.3,
(E T1 E T2) (zI z2) > (H(u ) H(u2)) (u1 u2),
(3.9)
1. 2
where u and u are solutions to SVIP[S(z ),H] and SVIP[S(z2),H],
1 2
respectively, and z1 7 z2. Since E has full column rank, u 1 u and
the strong monotonicity of H guarantees that the right hand side of
equation (3.9) is positive. Letting p = IIE1A 112 > 0, we have
(3.10) (H(u ) H(U2))(u u2)
2 a u u 2
2 2
= a lu u2
p
a
= 2
P
2
P
a
2
P
a
2
P
IIElA ul u2112
E1(Au1 Au2)112
1( z2
I(z. ) 2I2
1 2
where u E S(z ) and u2 S(z2).
obtain
(ET E T2 (zI 
Combine equations (3.9) and (3.10) to
z2~ 2 II z z2 I2'
P
3.3.1 A Special Case the TAP
Under the assumption that S is the feasible region of a TAP,
it is shown below that the mapping n(z) is point to point and Lipschitz
continuous.
34
For the TAP we have the following relationships:
S1 K 1 K k k k k k k
S' = { (w ,..., w y ,., y ) : Dw + B y = r ; x y w k 0 ,
k 1,..., K },
S 1 K k k k
S(w) = (x... x) B = Dw ; x 0, k = 1,..., K } where
1 K
w = (w ,...,w ), and
n(w) = { (DT ,..., DTK): C(X*) Bk = 0; A = 0, xk 0,
k = ,..., K } .
where, as before, C(F) = (C1(X),C2(Y)) is a nonnegative travel cost
vector.
For the kth OD pair, each element of the subvector D Tk is
the "time" to traverse the shortest path in the subnetwork which con
nects a pair of split nodes, i.e., nodes which appear in both the
master and the subnetwork. The following lemma shows that these times
are the same for all OD pairs.
Lemma 3.6 D k = D V k,A = 1,..., K.
Proof: From Theorem 3.1, Sk is the vector of dual variables for the
following minimum cost flow problem
min C1 (X)x
subject to
k k
If a and b are any pair of nodes in the subnetwork, then b a is the
b a
*
length of the shortest path from a to b, where [C1(X )]i, the ith
element of C1 (X ), is the cost associated with arc i. Also each
35
element of the vector DT has the form, ~ where a and b are
split nodes in the subnetwork.
Assume, without loss of generality, that for some particular
indices k and
[DT]i = a < [D ]i = E a.
Since [D Tk]i is the length of the shortest path from a to b in the
subnetwork,
[DT]i k a = k[C(X )]j,
kP
ab
where P is the set of indices for arcs belonging to a shortest path
ab
from a to b with respect to the kth OD pair. From the definition of
n(z),
[DTi] i a k([C(X*
[D i = = E ([C(X )]j A)
J ab
Sk Xj = E [C(X*)]j [D T ]i
j eP j eP
J ab J ab
= [D ]i [DT]i < 0.
This is a contradiction since X. > 0 V ij, 0
Theorem 3.7 If C (x) is strictly monotone, n(w) is a pointtopoint
mapping.
Proof: From Lemma 3.6, a vector in n(w) is completely determined by
any subvector D (T, for k e {1,...,K}. Furthermore, each component of
DT represents the shortest path time between pairs of connected split
36
nodes. Theorem 6.2 of Aashtiani and Magnanti [1981] guarantees that
these path times are unique. Thus, the vector D Tk is unique, i.e.,
the set n(w) is a singleton. O
Theorem 3.8 If C (X) is strongly monotone, then
x X* 11 II I I w w
X1 X2 "2 I l 2 112 '
*
where X1 and X2 are aggregate equilibrium flows to the traffic assign
ment problem for the subnetwork with the trip vectors, wl, and w2,
i.e., SVIP[S(wl),C1] and SVIP[S(w2),C], and B is a positive constant.
Proof: (Theorem 4.1, Dafermos and Nagurney [1982].) 0
Corollary 3.9 Theorem 3.8 holds for any standard norm.
Proof: The result follows from the comparability of any standard
norm. O
Theorem 3.10 If C1(X) is strongly monotone and Lipschitz continuous
with respect to the 1norm, i.e.,
C1(X1) C (X2) I11 p 11 x x2 Il p > 0
Then, n(w) is also Lipschitz continuous with respect to 1norm.
1%4
Proof: By Theorem 3.7, n(w) is a vector function (a pointtopoint
mapping). Define
*
(i) X. to be the aggregate equilibrium flow vector with
1
respect to the subnetwork trip vector, w., and
(ii) P to be a collection of arcs belonging to a utilized
pq
(shortest) path from p to q for the problem SVIP[S(wi),C ].
i'l
37
Then,
Iln(wl) n(w2 )1 = K
(p,q)
I [C (X1)]i 
i e
Pq
2[C1(X2)i I
i P
pq
where the outer summation is over the set of artificial arcs, i.e.,
pairs of connected split nodes, in the master network. Let
E 1[C (X1)]i 
iEPa
ab
S2[C(X2)]i
iEP
ab
=max { I [C (X1)]i
(pq) i P 1
pq
 2[C (X2) I ,
i ep
pq
and let 4 be the number of artificial arcs in the master network.
Therefore,
1In(w ) nl(w2) 1
S K E 1
ieP
[c1(X1)]i
 2[C (X2)]i
iE ab
ab
Without loss of generality assume that E 1[C (X1)]i
iEP .
. [C (X2 )]i,
ab
and
Ih(w1) n (w,) 1i 15 K.? {
 [C (X2)]j
jPab
ab
S1[C (X) ]i
i CP
ab
< K.* Z {
jeP
ab
5 K E 2
jeP
ab
[C1(X1)]j [Cl(X2)]j
S[C1(X1)]j
 [C1(X2)]j
KK p I II
5 K I p Xl X2 1
5 K p 8w w2 l1
38
In the above sequence of inequalities, the second inequality holds
because P2 is not necessarily the shortest path between a and b with
ab
respect to the cost vector C (X ); the third inequality holds because
2
for all i c P
ab
*
[C (X1)]i [C (X2)]i ; [C1(X1)]i [C1(X2)]i ,
the fourth inequality holds because
E 2 [C(X 1 )]j [C1(X2)]j  [C (X1)]j [C1(X2)lj
ePab jL1
= [C(x) Cl(X2) II
1 II
I1 2
and the last inequality follows from Theorem 3.8. Since all scalars
K, 0, p and 8 are positive, the theorem is proved. O
Corollary 3.11 Theorem 3.10 holds for any standard norm.
Proof: The result follows from the comparability property of any
standard norm. O
3.4 Algorithms for the Master VIP
Most algorithms for VIPs require the mapping to be at least
a continuous pointtopoint mapping. Thus, the discussion in this sec
tion will be limited to cases where the mapping n(z) is a vector func
tion, i.e., the traffic assignment problem.
1 K
It is clear that a feasible flow vector (x ,..., x ) with
respect to a subnetwork trip vector cannot be feasible to the subnet
work with respect to another subnetwork trip vector. Thus, the sets
S(") and 2 b 2
S(w ) and S(w ) would be disjoint when w r w Also, assume that
39
C1(X) and C2(X) are strongly monotone, C1(X) is a Lipschitz continuous
vector function, and C (Y) is a continuous vector function. By Corol
lary 3.4, Theorems 3.7 and 3.10, the mapping n(w) is a Lipschitz con
tinuous and strictly monotone vector function, and the composed mapping
(n(w),C2(Y)) is a continuous and strictly monotone vector function. In
this case, the cutting plane algorithm [Zuhovicki et al. 1969 and
Nguyen and Dupuis 1981] can be applied to the master problem.
Algorithm
'kO 0
Step 0 Let (w ,Y ) E S' and set t = 0 .
Ai i
Step 1 Define Pt(w,Y) = min (n(w (w w) + C2(Y) (Y Y)
i = 1,..., t) and let
t+l t+l
P3.8 Pt(w Y ) = max Pt(w,Y) : (w,Y) S' }
Step 2 If (w Y ) = (w Yj) for some j = 1,..., t, stop and
".t+l t+l
(w Y ) is the solution. Otherwise, set t = t+l and go
to Step 1. 0
In Step 1, the value of n(w ) is obtained by solving the sub
variational inequality problem for each w and problem P3.8 can be
reformulated as a linear programming problem [Nguyen and Dupuis 1981]
max Z
subject to
ZO < nT(w)( w w) + C2(Y)(Y1 Y) i = 1,..., t
(w,Y) E S',
where S' is the set of linear constraints described earlier.
40
Depending on the assumptions imposed on the original problem,
other algorithms mentioned in Chapter 2 can be applied to the master
problem as well. For example, if the subnetwork matrix, B1, is a node
arc incidence matrix of a tree, then n(w) is strongly monotone whenever
C1(x) is strongly monotone. With the tree structure, the matrices B
and D have linearly independent columns, and one of the rows in each
matrix is redundant. By deleting the same row from both matrices, S(w)
is equivalent to
1 K k ^k k
{ (x..., x) : x= Dw x 0, k 1,..., K.
where B1 and D are matrices B1 and D with one row missing, and they are
both nonsingular. Thus, the assumptions of Corollary 3.6 are fulfilled
and the above claim is proven. In the case of a strongly monotone
travel cost vector, most algorithms in Chapter 2 would be applicable to
the master problem whenever the subnetwork is a tree.
3.4.1 The Standard Traffic Assignment Problem
When C(F) = (C1(X),C2(X)) has a diagonal and positive semidef
inite Jacobian, the master and subvariational inequalities reduce to
the following pair of optimization problems of Barton and Hearn [1979]:
P3.9: (Master problem)
min V(w) + E Z.(Y.)
jcL2
subject to
k k k
Dw + B2y = r k = 1,..., K
k k
y ,w k 0
k
Y= Ey
41
P3.10 (Subproblem)/
V(w) = min E Z.I(Y
jcL 3
jEL1
subject to
k k
Bx = Dw k = 1,..., K
X = x k
k
x" 0,
where z.(X ) = / C (t)]. dt and Z.(Y) = Yj [C2(t)] dt. Barton and
I 1 0 '1 j J ) 0 2
Hearn solved both problems by the FrankWolfe algorithm and noted that
the time consuming step is the line search for the master problem.
Each evaluation of V(w) requires solving the subproblem, P3.10, which
is another standard traffic assignment problem. To alleviate this
problem, they suggested a heuristic method that performs the line
search, not on the objective function of the master problem, but on a
function defined to be the maximum of all tangent planes generated in
prior iterations. Below, we describe a convergent algorithm based on
the suggestion given in Geoffrion [1972].
From duality theory,
Sk k k
V(w) = maK i f { E Z (X ) u (B x Dw ) } ,
u x ?0 jEL1 k
where u is the Lagrange multiplier associated with the flow conserva
tion constraints for the kth OD pair. Equivalently, the master problem
can be stated as follows [Geoffrion 1972]:
min Z
subject to
SE Z.(Y.) + if { E Z.(X) uk B x Dw) V uk
jeL2 x 0 j EL k
2
42
k k k
k k k
y ,w 0; Y = E y,
and in practice the following relaxed version is used:
min Z
subject to
Z E Z.(Y.) + if { E Z.(X ) uk(B1xk Dw }
j0 EL2 x3 jeL 3 k
t= 1,..., T
k k k
Dw + B2y = r k = ,..., K
kk k
y ,w k 0; Y = E y
where ukt is usually the optimal Lagrange multiplier with respect to
the tth subproblem.
Algorithm [Geoffrion 1972]
Step 0 Let (w ,Y0) be an element of S'. Solve the subproblem V(w)
to obtain the Lagrange multiplier, uk'0 for all k. Set T = 0
0 '"0
and OBJ = E Z (Y.) + V(w ).
jEL2 Y 3
Step 1 Solve the following problem:
P3.11: min Z
subject to
k,t k
Z 2 E Z.(Y.) + E u (Dw
jEL2 k
+ if { E Z (X.) ukt(B x) } t = 0,..., T
x 20 jEL
43
k k k
Dw + B2y = r
k k
y ,w 0.
rV* *
Let (ZO,w ,Y ) be an optimal solution. If OBJ Z0 < ,
stop.
4 4* k,T+1
Step 2 Solve the subproblem V(w ) and obtain u .Set T = T + 1
and go to Step 1. 0
The property below makes problem P3.11 more appealing.
k,t 1 ""t k,t
Property 3.12 If xt, k = ,..., K solve V(w ), then x V k
also solve
P3.12: inf { E Z.(X ) I uk (B x ) }
xK 0 jeL1 k
kt
Proof: The vector uk must satisfy the following KuhnTucker condi
%It
tions for the subproblem V(w ),
C1( E xt) + Btu xt = 0 k = 1,..., K
k
Bk,t = Dk,t k = 1,..., K
k,t k,t) =
(X )(x )=0
Xk't,xk't Z 0
However, the KuhnTucker conditions for problem P3.12 consist of only
k't
the first and the last two sets of equations. Therefore, xkt must
also be a KuhnTucker point for problem P3.12, and by sufficiency of
the conditions xkt is also optimal. O
Define
tk,t k k,t k,t
6 = if { Z Z.(X ) E u (Bx ) } = E Z(X ) u Blx ,
x 0 jeL1 3 k jeL1 k
44
kt k t k t" t k,t
and by Lemma 3.6, E uktDw = uD E w = u DW where uD = u D for
k k
some k. Thus, the first set of constraints in P3.11 becomes
Z E Z.(Y.) + utDW + t t = 0,..., T
jeL2 3 3
and P3.11 is equivalent to
P3.13: min VT (W) + E Z.(Y.)
J L2
subject to
k k k
Dw + B2Y = r k = ,..., K
k k
.W = E w ; Y = E y
k k
k k
w ,y 0O,
where V (W) = max { utDW + 8t : t = 0,..., T } and P3.13 could be
solved by any method for the standard TAP.
To illustrate, consider the symmetric traffic assignment
problem in Figure 32 [Barton and Hearn 1979]. Figure 33 gives two
possibilities for decomposing the original network. The resulting
subnetwork in Figure 33(a) is a tree, and since there is only one path
between any pair of nodes, the solution to the subproblem is trivial.
As for the subnetwork in Figure 33(b), the subproblem is a symmetric
traffic assignment problem, and in our implementation we "solve" it by
performing five FrankWolfe iterations. For both choices of the sub
network, the approximate (or minmax) master problem, p3.13, is
"solved" by performing ten FrankWolfe iterations. Table 31 summa
rizes the results of the above computations. For comparison, it takes
the FrankWolfe algorithm fifty iterations or 1.84 CPU seconds (on the
IBM 3033N computer at the University of Florida) to achieve the
optimal flow
/
3 4
1 10 20
2 30 40
Trip Matrix
Zi(Fi) = F Ti(1+0.15(t/Ci) )dt.
Cost Function
Cost Function
Arc Ti Ci
1,5 5 10
1,6 6 16
2,5 3 35
2,6 9 .18
5,6 1 50
5,7 5 25
5,9 2 35
6,5 1 50
6,8 5 25
6,9 2 35
7,3 3 25
7,4 6 24
7,8 1 50
8,3 8 39
8,4 6 43
8,7 1 50
9,7 2 35
9,8 2 25
4
Volume Delay = T. ( +0.15(F./C.) )
Optimal objective value 1453.
Figure 32: An Example
45
23.63
46
(b)
Figure 33: Choices of the Subnetwork
47
Table 31: Computational Examples for the Minmax Master Problem
No. of Iter.
Figure 32(a)
Obj.
2950.00
1824.44
1556.91
1480.21
1463.87
1456.66
1456.22
1455.85
1455.70
1455.60
1455.55
1455.10
1454.97
Time*
0.17
0.50
0.83
0.97
1.30
1.37
1.44
1.64
1.81
1.88
2.56
2.99
Figure 32(b)
Obj .
2950.00
1465.51
1461.41
1459.98
1459.31
1459.13
1459.11
1459.10
1459.09
1459.00
1458.88
1458.82
1458.70
Time*
0.37
0.65
0.93
1.17
1.29
1.33
1.38
1.47
1.52
1.59
1.89
2.14
* The times are in CPU seconds on the IBM 3033N.
48
objective function value of 1457.47. However, one would not employ
transfer decomposition with the objective of improving computational
times the methods of Chapter 4 are more efficient. The purpose of
this example is to indicate that the technique is implementable when
necessary.
3.5 Applications: Link Abstraction and SubArea Focusing
Both link abstraction and subarea focusing are concerned with
replacing the original network with a simplified (and smaller) network.
Link abstraction refers to a procedure of collapsing a set of detailed
links into another (smaller) set of aggregate links. The main goal is
to abstract from the physical street or road network a simplified, but
accurate, representation to be used for analysis. In subarea focus
ing, the transportation planner is interested in a set of proposed
changes to street network within a given area. The accurate, but expen
sive, approach is to solve the entire network problem for each of the
proposed changes. The less expensive, inaccurate, and often used ap
proach is to discard the arcs external to the given area, estimate a
subnetwork trip table, and solve the traffic assignment subproblem for
each proposed change. A compromise is to replace the external arcs
with a simplified set of aggregate arcs.
Thus, both problems deal with the same type of simplification of
the original network, the physical street network for link abstraction
and the external arcs for the subarea focusing problem. The major
obstacle in both problems involves the process of incorporating the
original link attributes into the aggregate links. Many heuristic
methods [Chan 1976 and Hearn 1978] to construct the aggregate cost
49
functions for the aggregate links have been proposed, but none can
guarantee that the obtained solution reasonably approximates the true
solution.
Under the assumption that the travel cost vector is monotone
and diagonal, transfer decomposition provides a formal mechanism for
incorporating the original link attributes into the aggregate links.
In this case, the master and the subproblem are optimization problems,
i.e., problems P3.9 and P3.10. However, the network defining the two
problems are different under this application. Consider the network in
SFigure 3.4(a) as the entire network, and let the indicated part of the
network represent the area to be kept unaltered; this might represent
the network of principal arterials for the link abstraction process or
the target area for the subarea focusing problem. We shall refer to
this unaltered area as the principal area of the network. It is clear
that the principal area must be kept in the master network since it
cannot be replaced by pseudoarcs (the set of aggregate arcs) and the
rest of the network is left in the subnetwork(s). Thus, the resulting
problems are
P3.12: (The master problem)
min Z.(X.) + V(w)
j EL3 3
subject to
k k k
B x + Dw = r k = 1,..., K
X = xk
k
k k
x ,w O0
k = ,..., K.
50
Principal Area
(a) The Entire Street Network
o

0'
\ /
S/ /
\ /
\ /
/ \
/ \
/
/ \
i, D
(b) The Simplified Problem
Figure 34: An Example of Network Simplification
51
P3.13: (The subproblem)
V(w) = min E Z.(Y.)
iL2 1 1
subject to
k k
B2y = Dw
Y= E yk
y 0
k = 1,..., K
where B1 and L are the nodearc incidence matrix and the arc set, re
spectively, for the principal area. Similarly, B2 and E2 are the node
arc incidence matrix and the arc set, respectively, for the network
k
external to the principal area. As before, w and D are the artificial
vector and its nodearc incidence matrix. Furthermore, Barton and
Hearn [1979] showed that the objective function V(w) is also a func
tion of the aggregate flow on artificial arcs; i.e.,
P3.14: V(w) = V(W) = min Z.(Y.)
subject to
B2y(p,q) = Dw(pq) V artificial arcs (p,q)
Y = Z y(p'q)
y(P'q) r 0 V artificial arcs (p,q),
k
where W = E w Note that the OD pair index is different from before.
In the subproblem, the origin and destination nodes are the split nodes
and any of the original origin and destination nodes already included
in the subnetwork.
Since V(W) depends on the aggregate flow on the artificial
arcs, it must be the exact aggregate cost function for the aggregate
52
links in the above discussion. If V(W) is known explicitly, solving
P3.12 with V(w) replaced by V(W) would produce the correct aggregate
A
flow in the principal area, i.e., X Unfortunately, V(W) is defined,
in terms of an optimization problem. There is, however, one very
A
special (and obvious) case where the function V(W) is explicitly known.
Specifically, assume that all external arcs have linear cost. Then, as
illustrated in Figure 35, linear path costs could easily be assigned
to the artificial arcs. These path costs are the length of the short
est paths which connect the head and tail nodes of the artificial arcs
and consist of only arcs.external to the principal area.
Since the evaluation of the function V(W) generally requires
solving an optimization problem, any algorithm which requires many
evaluations of V(W) would be quite slow. One way to alleviate this
problem is to approximate V(W). One of the two easily implementable
A
methods is curve fitting and the other is using the function V (W).
A
As defined earlier, V (W) is the maximum of tangent planes generated by
the algorithm in the previous section. For curve fitting, one can
obtain different approximations of V(W) by varying the form and the
norm used in calculating the error of the approximate function. For
the example problem below, we consider the following four variations of
the curve fitting scheme:
(i) LSQ1: V(W) = a + bW + WAW.
Norm = Euclidean norm.
(ii) LSQ2: V(W) = b + WAW.
Norm = Euclidean Norm.
(iii) MINMAX1: V(W) = a + bW + WAW.
Norm = Chebyshev norm.
53
Cost function
\
Principal Area
(a) The Original Network
Cost function
2W1,7
W9,3
'>
I/ 10U
 >(T)
^lO. v
Principal Area
(b) An Equivalent Network
Figure 35: A Network with Linear External Arcs
54
(iv) MINMAX2: V(W) = bW + WAW.
Norm = Chebyshev norm.
where a is a scalar, b is a vector and A is a matrix.
Consider the symmetric traffic assignment problem in Figure
36. Assume that the original network is decomposed as in Figure 37,
then the subproblem consists of two independent symmetric traffic
assignment problems, denoted as V (W) and V (W). For the four curve
fitting schemes, we randomly generated fifty flow vectors for the
artificial arcs from the set { (W1, W2, W8) : 0 W1, W2 5 30,
0 W W3 W4 70, 0 ~ W5,W7 S 40, 0 S W6, W8 : 60 } using a uniform
random number generator. For the second approximation method, the
A A
functions, V T(W) and V T(W), used in the calculations below are
obtained by solving the decomposed problem using the method described
in the previous section where T is equal to twenty, i.e., V i(W) is
taken to be the maximum of the first twenty tangent planes generated by
the algorithm.
Table 32 displays the resulting flow vectors for arcs inside
the principal area and the errors due to each approximation. From
these results, it is encouraging that the errors in the solution are
reasonably small. Also, it is interesting to note that V (W) achieves
the smallest errors among the five approximations.
In spite of the fact that the correct aggregate function is
generally unavailable in a convenient form, the analysis in this chap
ter provides the first rigorous basis for construction of an aggrega
tion function in practice. The construction of suitable approximate
functions along the lines suggested above is a matter for extensive
computational experimentation and is left to future research.
55
3 4
1 10 20
2 30 40
Trip Matrix
F 4
Zi (F) = Ti( +0.15(T/Ci) )dt
0
Cost Function
Optimal Objective Value = 2750.86
Figure 36: Original Network
Arc
1 5
1 6
2 5
2 6
5 6
5 7
5 8
6 5
6 7
6 8
7 8
7 9
7 13
8 7
8 10
8 13
9 10
9 11
9 12
10 9
10 11
10 12
11 3
11 4
11 12
12 3
12 4
12 11
13 9
13 10
56
S(w 
1r
(W)
1
 '
>
\w
Figure 37: The Decomposed System
A
V (W)
>
57
Table 32:
Computational Results
ARC OPT. FLOW LSQ1 LSQ2 MINMAX 1 MINMAX 2 V (W)
T
7,8 5.84 0.0 0.0 0.0 9.16 0.0
7,9 23.55 22.85 22.33 25.97 20.42 24.94
7,13 47.99 43.13 44.22 47.32 40.83 48.61
8,7 0.0 0.0 0.0 0.0 0.0 0.0
8,10 3.25 0.0 0.0 0.0 9.16 0.0
8,13 25.21 34.02 33.46 26.71 29.58 26.45
9,10 0.0 0.0 0.0 0.0 0.0 0.0
10,9 0.0 0.0 0.0 3.61 0.0 0.0
13,9 40.25 45.33 43.45 45.68 44.18 43.62
13,10 32.94 31.82 34.22 28.35 26.24 31.44
Abs. error in the soln.
(1) Euclidean norm
(2) Chebyshev norm
% error in the soln.
(1) Euclidean norm
(2) Chebyshev norm
Error in obj. function
(1) Absolute error
(2) Percent error
13.17
8.81
11.84
8.25
16.67% 14.99%
18.36% 17.19%
97.68
3.55%
254.8
9.26%
10.81
5.84
13.68%
12.16%
562.57
19.14%
13.65
9.16
17.28%
19.08%
267.40
9.72%
7.20
5.84
9.11%
12.16%
30.16
1.09%
CHAPTER FOUR
APPLICATIONS OF SIMPLICIAL DECOMPOSITION
Simplicial decomposition was originally developed for a
convex programming problem with linear constraints, i.e.,
P4.1: min { Z(x) : Bx = b, x L 0 }
where Z(x) is a convex function, B is an m x n matrix, b is a constant
vector in Rm, and x is a.vector in Rn. To simplify the notation, the
feasible region S = { x : Bx = b, x a 0 } is assumed to be bounded, and
thus it is compact and representable as a convex hull of its extreme
points. Although the number of extreme points of S is finite, but
generally large, making it impractical to enumerate them a priori.
Therefore it is more efficient to obtain a small subset of extreme
points initially and generate others as required. This is exactly the
basic scheme of simplicial decomposition which can be precisely stated
as follows:
The Basic Algorithm
0 0 ^ 0
Step 0 Let x be an element of S. Set S = {x } and x = x .
Step 1 Let y solve the subproblem: min { VZ(x)y : y S and
set n = n U {y).
Step 2 Let x solve the master problem: minimize Z(x) over the con
vex hull of 2. If x solves problem P4.1, stop. Otherwise,
purge 1 of all extreme points with zero weight in the expres
sion of x as a convex combination of extreme points in Q. Go
to Step 1.
58
59
The advantages of this basic algorithm are (i) all con
straints are assigned to the subproblem; thus the iterates are always
primal feasible and no dual variables are necessary to link the master
to the subproblem, (ii) the number of extreme points is kept minimal by
dropping those with zero weight, (iii) the master problem is relatively
simple, and (iv) its convergence rate coincides with that of the algo
rithm employed for the master problem. From a theoretical standpoint,
the disadvantages are the upper bound on the cardinality of 0 during
any iteration is n+l (Caratheodory's Theorem, e.g.; Bazaraa and Shetty
[1979]) and the master problem must be solved exactly. The latter does
not pose any real difficulty. In practice, the master problem can
never be solved exactly on a computer and the method still produces
good solutions. The requirement for the exact solution originates from
the argument often used to establish convergence. It rests on the fact
that the number of extreme points is finite and that to each simplex
one can associate with it the optimal objective function value of the
corresponding master problem. Since the optimal objective function
value decreases from one iteration to the next, the same simplex can
never be regenerated and the algorithm must terminate after a finite
number of iterations. However, there is an alternative argument
[Holloway 1974] which only requires a decrease in the objective func
tion value (not necessarily optimal) to insure convergence. Using this
argument, the algorithm proposed below only requires an approximate
solution to the master problem.
As for the first disadvantage, it is impractical, if not im
possible, in most applications to reserve enough computer memory to
store n + 1 extreme points. Often, one would allocate enough computer
60
memory for only a few extreme points and hope that it is sufficient.
Below, we propose a modified simplicial decomposition scheme, called
restricted simplicial decomposition, which allows the user to restrict
the number of extreme points to any desired number r Z 2. Computa
tional experiments demonstrate an interesting and potentially very
useful tradeoff between r and the progress of the algorithm.
After the treatment of restricted simplicial decomposition,
the same basic decomposition scheme is extended to the VIP. A straight
forward application of the basic algorithm to the VIP does not work; a
counter example showing that the resulting algorithm cycles through
three extreme points is given in Magnanti [in press]. Most authors
[Bertsekas and Gafni 1982 and Smith in press a] circumvent the problem
by disallowing column (extreme point) dropping. In the algorithm for
VIP below, column dropping is allowed when the algorithm has made sig
nificant progress toward the solution, and, as before, the exact solu
tion to the master VIP problem is not required.
Before proceeding, it is important to note that the set of
feasible flow vectors for the TAP is an unbounded polyhedral set ex
pressible in the form of the feasible region of problem P4.1. For a
TAP with a nonnegative travel cost function C(*), it is easily seen
that the only flow vectors of interest are those that lie in the convex
hull of all extreme aggregate flows. An "extreme aggregate flow vector
(or pattern)" is so called because it is an extreme point of the set of
feasible flow vectors and is the sum of individual flow patterns which
correspond to minimum cost paths between OD pairs for some given value
of C(*). Since the convex hull of these extreme aggregate flows is
compact, the simplicial decomposition scheme is applicable to the TAP.
61
Henceforth, the feasible region S for the TAP is taken to be the above
convex hull.
4.1 The Restricted Simplicial Algorithm
This algorithm is a variation on the basic algorithm present
ed earlier for problem P4.1. It may be viewed as an application of
the projection technique developed by Goldstein [1964] and Levitin and
Polyak [1966] to the master problem, combined with a scheme for retain
ing at most r grid points where 2 5 r n n + 1 is specified by the user.
Grid points are defined to be elements of S which define the simplices
in Step 2. In the basicalgorithm, grid points are taken to be extreme
points here they may be prior iterates (the x) as well. In fact, it
is crucial to the convergence of the method that the immediate prior
iterate be in the convex hull of the defining grid points in the next
iteration. Without this feature, monotonicity of the objective
function could not be guaranteed.
Although there are many methods for solving the master prob
lem, the projection technique is a natural choice. It allows for the
approximation of the master problem solution, and the accuracy can be
controlled by varying the number of projections performed per itera
tion. Also, the simplicity of the master problem's constraints permits
an almost closed form solution to the projection problem.
Bertsekas [1979] described the projection algorithm as
follows: assume that Z(x) in problem P4.1 is a twice continuously dif
ferentiable convex function. (For the TAP, the function Z(x) is as de
fined in equation (2.17) and when C(*) is monotone and continuously
differentiable the function Z(x) is convex and twice continuously
differentiable.) Since the set S is assumed to be compact, there
62
exists a constant Y > 0 such that
(4.1) IlVZ(x) VZ(y)  ; y i x yl V x,y E S,
i.e., VZ(x) is also Lipschitz continuous [Bertsekas 1976]. In addi
tion, assume that {Mk) is a sequence of matrices satisfying
(4.2) a llx 112 < x Mkx s5 lx 12 V x null space of B, and V k
where a and 8 are positive scalars and
(4.3) a > y/2.
Then an approximation of the objective function using Mk can be
written as
(4.4) Z(y I ) = Z(x) + VZ(x)(y x) + ((y )Mk(y x)
and the solution to
P4.2: min { Z(y x) : y S
is the projection of the vector x Mk VZ(X) onto the set S with
respect to the norm lIx I= (x Mx) Under the given assumption, Bert
A
sekas [1976] proved that x's generated successively by solving problem
P4.2 converge to x the solution of problem P4.1.
Since problem P4.2 can be formulated in terms of extreme (or
grid) points of S, i.e.
P4.3: min { Z(AX x) : 1. = 1, X. > 0 V i }
1 1
where columns of A represent the extreme (or grid) points of S, its
solution can be calculated quite easily when the matrix A MkA is
k
63
diagonal. The solution method is almost closed form because the objec
tive function of problem P4.3 is a strictly convex quadratic in the A
variables with no crossproduct terms. The required calculations are
similar to those described in Section 5 of Held et al. [1974].
Let 0 denote the set of retained grid points and H(0) denote
the convex hull of f, i.e.,
H() = x : x = Z X.a E A. = 1, X. i 0 and a e }
Restricted Simplicial Decomposition Algorithm
Step 0 Let x be an element of S. Set 0 = {x0} and t = 0.
Step 1 Solve
P4.4: VZ(xt)yt = min { VZ(xt )y : y S }.
If VZ(x )(y x ) 0, stop and x is a solution. Other
wise,
(i) let St+l = t U {yt) if It I < r, or
(ii) let nt+l be St with two of its elements replaced by x
and yt if I t = r.
L
Step 2 Solve
P4.5: Z(xt+lI t) = min { Z(x t) x E H(t+l) }
Purge t of all grid points with zero weight in the expres
t+l
sion of x as a convex combination of the grid points. Set
t = t+l and go to Step 1. O
It is important to observe that x always belongs to H(t+l)
at the beginning of Step 2. In Step 1, if it I < r, then xt H(Qt)
implies that xt H( t U {yt} ); otherwise, x belongs to
64
H( t+) by construction. This observation is later used to prove that
the objective function decreases at each iteration.
Under the column dropping scheme in Step 1, the set of re
tained grid points always consists of the newly generated extreme point
and the previous iterate for the case r = 2. This is exactly the
scheme in Frank and Wolfe [1959]. When r > 2, two of the retained grid
points are replaced by the same two points mentioned above if the total
number of retained grid points equals r. Otherwise, the newly generat
ed extreme point is added to the current set of grid points.
We now prove the convergence of the restricted simplicial
algorithm. The argument is a synthesis of known results from Bertsekas
[1978] and Holloway [1974]. Lemmas 4.2 and 4.3 show that the approxi
mate objective function decreases from iteration to iteration which, in
turn, implies that the same is true for the real objective function.
From this conclusion, the convergence proof follows.
*
Lemma 4.1 VZ(x ) (y x ) = 0 V y e S if and only if x is an opti
mal solution to problem P4.1.
Proof: This is a standard nonlinear programming result. O
Lemma 4.2 If x is not optimal to problem P4.1, then
Z(xt+l x) < Z(x x = Z(xt).
t t+l t+l
Proof: By the above discussion, xt H( t). Since x minimizes
t t+l
Z(x I xt) over H(tl)
t+1 t t t
Z(x t xt (x xs
If Zt+l It t
If Z(x Ix ) = Z(x x ), then x is also a minimum and
65
VZ(xt [ x) (y x) = z(xt) (y xt) 0
which contradicts the fact that x is not optimal to problem P4.1 and
the stopping criterion is not satisfied, i.e., VZ(xt) (yt x ) < 0. O
Lemma 4.3 Given that Mt, t = 1, 2,..., satisfy the assumptions above,
t+1 I t ^
Z(x I x) < Z(xt xt) implies Z(xt+) < Z(xt)
Proof: The argument is essentially from Bertsekas [1979]. Since
t+l
x solves problem P4.5,
St+1+ t t+t t
VZ(xt+l xt) (x xt+l) 0 V x H(t)
(4.5) or [ VZ(xt) + Mt(xt+ x ) ] (x x tl) a 0 V x e H() .
In particular, let x = x and (4.5) becomes
t t+1 t t+1 t
SVZ(x ) + Mt(xt x ) ] (x xt ) 0.
(4.6) VZ(xt)(xtl x ) (x x )Mt(x x )
The second order Taylor expansion of the objective function at the
t
point x gives
Z(xt+) Z(x ) + VZ(x ) (xt+l x)
1
+ / (VZ(x + u(x x)) VZ(xt))(xt+l x)du
0
t t+1 t t+1 t
Z(xt) Z(xt+l) = VZ(xt)(x+ xt)
1
t t+l t t t+1 t
f (VZ(xt + u(x x)) VZ(x ))(x tl x )du
0
66
t t+l t
SVZ(xt) (x t+ x )
II Z(x + u(x t+ x ) Z(xt) II Ix xt II du
0
1
S(x t )Mxt+1 x f u Y x +l t 2 du
0
Sa IIxt xt2 /2 x xt 112 = (a /2)Ixt+l xt2 > 0.
The first inequality follows from the Schwartz inequality, the second
follows from the Lipschitz condition and (4.6), the third follows from
t+l t
the property of Mt and the last follows from the fact that x + x O
Having established the monotonicity in the objective function
values of the iterates, the convergence of the sequence {xt) to an op
timal solution of problem P4.1 follows from arguments similar to those
of Lemma 4.2 and Theorem 4.3 of Holloway [1974] (see also Zangwill
[1969]). Therefore, the proofs for the next two results are omitted.
Lemma 4.4 A sequence {xt}t generated by the algorithm cannot admit a
subsequence {xt}t T satisfying the following properties:
(i) x x t T
(ii) y + y" t e T
(iii) VZ(xO")(y x) < 0.
Theorem 4.5 Under the stated assumptions, the algorithm converges to
an optimal solution of problem P4.1.
The algorithm as stated only performs one projection per
t 1 t
iteration, and since the projection of the vector x M Vz(x ) onto
H(t+ ) is not necessarily optimal to the original master problem,
i.e.,
67
P4.6 min { Z(x) : x e H(St+l) }
the same extreme point, y is often successively regenerated in Step
1. One way to insure that yt does not belong to t is to solve the
original problem exactly. However, the following theorem suggests that
only an approximate solution is required.
Theorem 4.6 Let x be a solution to problem P4.6, but not a solution
to problem P4.1, and 4 be the set of solutions to min { VZ(x)y : y S }.
Then there exist a 6 > 0 and a yE # such that y solves
min { VZ(x)y : y S I and Vf(x)(y x) < 0 for all x satisfying
lix x '5 6.
Proof: Since VZ(x) is continuous, we have that for all e > 0 there
exists a 61 > 0 such that IIVZ(x) VZ(x) I< e whenever ix x < 61.
Furthermore, the problem min { VZ(x)y : y c S } is a linear program and
*
there must exist an e > 0 such that there is an element y from the set
I which solves min { VZ(x)y : y E S I whenever IVZ(x) VZ(x) < e ,
or whenever x x I < 61'
Define
g(x) = VZ(x)(y x),
then g(x) is continuous and g(x) < 0, since x is not a solution. Hence,
*
there exists a 62 > 0 such that g(x) < 0 for all x satisfying
Ix x I 6 To obtain the result, let 6 = min{ 1 6 .
2 1' 2
Therefore, preventing the regeneration of the same extreme
point requires a good approximate solution to the master problem which
can be obtained by performing many projections per iteration. That is,
Step 2 can be replaced by the following steps:
68
0 t
Step 2a Let w = x k = 0 .
Step 2b Solve
k+1 k w k t+1
Z(wk+ w ) = min { Z(w w) : H(l) }
k+l
Step 2c If w does not satisfy a stopping criterion, set k = k + 1
t+l k+l
and go to Step 2b. Otherwise, set x = w and purge
t+l of all grid points with zero weight in the expression of
t+i
x as a convex combination of grid points. Set t = t + 1
and go to Step 1. O
k.
As long as w is not optimal to P4.6, Lemmas 4.2 and 4.3
guarantee that
k+1 k
Z(wk+l) < Z(w ) V k = 1, 2,...,
which further implies that
Z(xt+1) < Z(x .
From this inequality, the results of Lemma 4.4 and Theorem 4.5 follow.
4.1.1 Computational Testing
Since the FrankWolfe direction finding problem, i.e., the
problem in Step 1, is generally nontrivial and the projection is rela
tively easy to perform when A MkA is diagonal, the following stopping
rule of Step 2c is adopted:
%X t+l I t+1
Stop the projection process if G(w I ) 5 e or
t+l > MAXP, where
SII ut+1 '
(i) G(x ) = VZ(x)(x y)/ VZ(x)y, and y solves
min { VZ(x)y : y e H(k+l) .
69
(ii) MAXP = maximum allowable number of projections.
The numerator of G(x I t+l) is the gap function [Hearn 1982] at x with
t+l I t+l
respect to the region H(n ), and G(x 0l ) is therefore the gap
function normalized by VZ(x)y. We refer to it as the relative gap
function. Instead of letting et be constant for all t, the sequence t
is defined as follows:
e = min { e ; 6* G(xt S) }
st tmin{ t_l;1
where 0 < 6 < 1 and EO is set a priori. This construction allows the
values of t to be relatively large in the early stages of the algo
rithm since the main interest at this stage is to generate grid points
and not to solve the master problem. However, in the latter part of
the algorithmic process, the algorithm presumably would have already
obtained the simplex which contains a solution, and a good approximate
solution to the master problem would mean a good approximate solution
to problem P4.1. So, by making the et small at the latter part of the
algorithm, more effort is put into solving the master problem. The
constant MAXP is only used to prevent an excessive number of projec
tions during each iteration.
The test problems include three traffic assignment problems,
the symmetric problem (Figure 41) from Nguyen and Dupuis [1981], the
nine node problem (Figure 42) and the windowpane problem (Figure 43)
from Barton and Hearn [1979]. For the experiments the matrix AT MA is
taken to be the diagonal matrix whose diagonal elements all equal to
k
.75 times the gap value at w .
The program for restricted simplicial decomposition is writ
ten in FORTRAN with variables r, 6 and MAXP as its parameters. Ini
tially, a (free) PDP 11/34 computer was used to run the program for
70
3 4
1 400 800
2 600 200
Trip Table
Figure 41: Nguyen & Dupuis Symmetric Problem
71
S1 10 20
9
S/ 2 30 40
6 < 8 4
Trip Table
Figure 42: Nine Node Problem
8 4
1 2 3 4
1 10 20 10
2 10 20 10
3 10 10
4 10 10 20
Trip Table
Figure 43: Windowpane Problem
72
each problem and each r (the number of grid points) at various values
of the parameters 6 and MAXP in order to determine their "optimal" val
ues, 6 and MAXP For each pair of 6 and MAXP values, the program was
run until it achieved an objective function value which is as good as
or better than the objective function value at the 100th iteration of
the FrankWolfe algorithm (written by Nguyen and James [1975]). For
each run, the following data were recorded: the number of iterations,
the number of projections performed, the CPU time and the final objec
tive function value. As an example, Table 41 summarizes the results
for the ninenode problem with r = 4. The CPU times in this table are
estimated CPU times required to execute the same program on the IBM
3033N computer (the PDP 11/34 does not have a CPU clock). Using these
estimated times as the sole criterion, Table 42 lists the best values
of 6 and MAXP for each value of r, the number of retained grid points,
for all three problems.
For the final comparison, we ran the programs using the para
meters listed in Table 42 on the IBM 3033N computer and plotted the
objective function values against the actual CPU times on the IBM
3033N. The results are summarized in Figures 44, 45 and 46.
From the results above it is clear that as r increases the
time required to achieve the objective function at the 100th iteration
of the FrankWolfe algorithm decreases. Furthermore, it is encouraging
that the value of r need not be too large to obtain fast convergence.
In all three problems, the largest value of r is 7, and this is quite
reasonable when compared to the Caratheodory bound 20 for the symme
tric problem, 19 for the nine node problem and 25 for the windowpane
problem.
73
Table 41: Computational results for the nine node problem with r = 4.
MAXP 100
6X
.1
.2
.3
.4
.5
.6
.7
obj. fun.
(I iter., I proj.)
Est. CPU time (sec.)
The objective function value at the 100th iteration of FrankWolfe algorithm
is 1455.9588 (3.10 sec.) and the best lower bound is 1453.1054.
1455.8660 1455.7278 1455.5480
(13,1025) (11,1091) (11,1252)
.6966 .6710 .7297
1455.8140 1455.7813 1455.8455
(16.1087) (13,1092) (12,1068)
.7937 .7210 .6874
1455.7771 1455.9149 1455.9149
(20,1124) (19,1047) (19,1047)
.9065 .8536 .8536
1455.8676 1455.9387 1455.9387
(31,1075) (30,1037) (30,1037)
1.1619 1.1232 1.1232
1455.9497 1455.9497 1455.9497
(42,996) (42,996) (42,996)
1.4053 1.4063 1.4063
1455.9374 1455.8856 1455.8856
(45,997) (43,1029) (43,1029)
1.4812 1.4432 1.4432
1455.9219 1455.9219 1455.9219
(67,1031) (67,1031) (67,1031)
2.0400 2.040 2.040
74
Table 42:
"Optimal" values of 6 and MAXP.
Nguyen & Dupuis Ninenode Problem Windowpane Prob.
Symmetric Problem
100th iter. of FU 100th Iter. of FW = 100th Iter. of FW 
85100.07 (3.12 sec.) 1455.91 (3.10 sec.) 1945.99 (3.41 sec.)
85098.5 1455.85 1945.98
(0.1,200) (0.1, 200) (0.1, 100)
(40,3804, 1.94) (29, 1699, 12.9) (47, 3898, 2.77)
85099.68 1455.72 1945.86
(0.1, 100) (0.1, 200) (0.1, 200)
(22, 1963, 1.28) (11, 1091, 0.71) (6 5T,07;'0.44
85097.5 1455.84 1945.7898
(0.2, 300) (0.1, 100) (0.4, 100)
(6, 223, 0.32) (9, 704, 0.55) (3, 44, 0.16)
85029.9375 1455.44
(0.1, 100) (0.1, 100) same as above
(5, 95, 0.22) (6, 171, 0.26)
1455.19
same as above (0.5, 100) same as above
(6. 118, 0.24)
objective function
(6, MAXP)
(I iter., # proj., time in sec.)
75
obj. fun.
87,000
86,50
86,000
85,500 r2 (Frank & Wolfe)
r3
r4
r5
85,000 r6,7
0.5 1.0 1.5 sec.
Figure 44: Nguyen & Dupuis Symmetric Problem
76
obj. fun.
1490. 
1480
1470
1460
2 (Frank & Wolfe)
1.5 sec.
Figure 45: Nine Node Problem
77
obj. fun.
1980
1970
1960 r=2 (Frank & Wolfe)
r3
1950
05 sec.
0.5 1.0 1.5 sec.
Figure 46: Windowpane Problem
78
In summary, the tests show empirically that there is a trade
off between computer storage used and the progress of the algorithm.
The implication for practical applications is that choosing r > 2 can
offer an improvement over the FrankWolfe method (where r = 2).
4.2 Simplicial Decomposition for the Asymmetric Problem
A simple extension of the basic algorithm to the VIP[C,S]
would be to replace the optimization problem in Step 2 with the corres
ponding VIP. However, as Magnanti [in press] has recently noted with
an example due to J. Hammond, the resulting algorithm cycles. To
prevent cycling, we onlyallow column dropping when the algorithm has
made satisfactory progress toward the solution. For the symmetric
case, the objective function was used to monitor the progress of the
algorithm. To establish the convergence of the simplicial technique
for the VIP[C,S], some choice of monitoring function must be made since
no objective function exists. There are several choices: the gap func
tion (G(F), Section 2.1), its dual function (g(F), Section 2.1), and
the user objective function [Smith 1983]. For simplicial decomposi
tion, the gap function is the natural choice since it is the byproduct
of solving the subprogram in Step 1. The gap value is simply the quan
tity C(x )x minus the optimal objective function of the subproblem.
In addition, it is as yet uncertain how the other two functions might
be evaluated computationally.
The algorithm of the VIP requires a small positive tolerance
6 and a convergent monotone sequence {t I where t > e > 0, and
St 0 as t t
et + 0 as t + co
79
The SDVI Algorithm
Step 0: Let x0 be an element of S, and define 0 = {x0 DO = 0 and
O
G = Set t = 0.
Step 1: Solve
C(xt)yt = min { C(xt)y : y e S .
If G(x ) = C(x ) (x y ) = 0, stop and x is the solution.
Otherwise,
(i) If G(xt) t t 6, let tl = t U yt
(ii) If G(x ) < t 6, let nt+l = (st\D ) U y
t+l t t
Let G = min { G G(x ) .
Step 2: Find xt+1 H(t+1) such that
t+l t+l t+l
'C(x )(x x ) 2 V x H(SI ) .
Let Dt+l be columns of 2t with zero weight in the expression
t t
of x as a convex combination of extreme points in t.
Increment t and go to Step 1. D
The parameters et allow for approximate solving of the master
variational inequality problem in Step 1. In fact, when C(x) is strong
ly monotone, i.e., there exists an a > 0 such that
(C(x1) C(x2)) (x x2) 2 a 1 x211 V x1, x2 S
At
and if x solves VIP[H( t),C] exactly, then
a lxt x2 < (C(xt) C(x))(xt x)
=C( t) (x x) C(x)(xt x)
80
t t ^.
5 C(x) (xt x)
Smax C(x) (x w) et
weH( t)
.lxt ^112 c/
Therefore, depending on a, the sequence et can be chosen to
t ^
assure x as near x as desired.
4.2.1 Convergence Proof
For the convergence proof we say that the algorithm completes
one major iteration every time the gap, G(xt), decreases, and we refer
to the iterations between any two major iterations as minor iterations.
Define T to be
t t t t
T = { t : C(x) (x y ) < G t = 1,2,3, ... } ,
then T is the index set of the major iterations.
It is clear that if the algorithm terminates at the kth
iteration, x must be the solution of VIP[S,C]. Also k e T since
C(xk (x y) = 0 < min { C(xt (xt yt) : t < k } and C(xt) (xt y
> 0 V t < k. Therefore, in the case of finite termination, the algo
rithm yields the solution at a major iteration. Otherwise, the algo
rithm generates a sequence {x }, and we will show that every convergent
subsequence of {x t ^ converges to the solution of VIP[S,C]. It
should be noted that the algorithm does not drop columns at every major
iteration; instead, it only drops columns when there is a sufficient
decrease in the gap, i.e., by an amount 6. As we will show, the para
meter 6 assures that ultimately no columns will be dropped, and this
leads to the proof of the following convergence theorem.
Theorem 4.7 Assume that the algorithm generates an infinite sequence
81
t
{x } and T is an infinite set, then any convergent subsequence of T
converges to the solution of VIP[S,C].
Proof: Since xt S for all t e T where S is compact, there must exist
t ^ *
a convergent subsequence {x } where T T. Let x be the limit
tCT
t
point of the subsequence {x t} Since the gap function, G(x), is a
nonnegative continuous function on S, there is an integer K such that
for all t e T and t > K
IG(xt) G(x*) < 6.
By monotonicity of the sequence {G(xt) tT' we have
t *
(i) G(x) G(x ) < 6 V t T and t > K
(ii) For any tl, t2 T such that t2 > t > K
t t2 t1 *
G(x ) G(x )
Note that (ii) also implies that for any t such that t e T and tl < t < t2
t t t
(iii) G(x 1) G(x ) < Gx ) G(x 2) < 6.
From (ii), (iii), and the fact that there is no column dropping during
any iteration t T there must be no column dropping after iteration K,
t
and the sequence of sets n t > K is a monotonically increasing se
quence. Since S has a finite number of extreme points, the sequence
t ^ *
{ }t>K must converge to some set 0, i.e., there exists an integer K
such that
t ^
0 = V t > K.
For all t,
C(xt) (x x) t V x e H( t)
82
ort t
max C(xt) (x w) 5 e
weH(t)
t ^ ^ t
However, since t = Q for t > K, y must be an element of
A t
^ = t where
C(x) (x y ) = max C(x) (x y).
yeS
Thus; for all t e K
t t t t t t
G(xt) = C(xt) (x y ) = max C(x) (x w) S .
wEH(Qt )
Because t is a sequence converging to zero,
t *
lim G(x ) = 0 = G(x ) .
tET
Therefore, x must be the solution of VIP[S,C]. Since the
subsequence defined by T is arbitrary, the proof is complete. 0
In the proof above we assume that T is an infinite index set
whenever the algorithm generates an infinite sequence. Now we show
that this assumption is in fact a true statement.
Lemma 4.8: If the algorithm generates an infinite sequence, then T is
an infinite set.
Proof: Assume otherwise; that is, there exists a K such that for all
t > K
t = C(xt t t k
G(x) = C(x) (x y ) > G = 8 > 0.
But, this implies that there is no column dropping after iteration K
and using a similar argument as in the above proof it can be shown that
G(x ) < E V t > K.
Since converges to zero, we obtain a contradiction.
Since s converges to zero, we obtain a contradiction. O
t
83
4.2.2 Discussion of Algorithm
In principle, any convergent method for solving variational
inequalities could be applied to the master problem of Step 2. How
ever, especially powerful methods, owing to the simplicity of the con
straint structure, are projection methods. Furthermore, Bertsekas and
Gafni [1982] have recently developed a convergent projection method
which is applicable to our master problem. Specifically, they have
shown that their method is linearly convergent whenever C(x) is Lip
schitz continuous and strongly monotone, even if there is a linear
transformation of the variables of the form x = AX. In other words,
the linear convergence result will hold if the problem is solved in the
space of X rather than of x. Their proof is important because the as
sumed properties of C may not carry over to the variational inequality
problem in the X space (see the discussion in Bertsekas and Gafni
[1982]). In all of our experimentation we have solved the master prob
lem of Step 2 in X space and defined Dt by those columns of ft for
which the associated X variable is zero.
Bertsekas and Gafni have also employed their method for sol
ving the traffic assignment problem. However, there are important dif
ferences between the strategy of their implementation and ours. Their
method works in the space of path flows, ours in the space of aggregate
extreme flows (as does the FrankWolfe method). Although, there are
generally more extreme aggregate flows than path flows, this does not
necessarily imply that more will be required in the expression of a
given flow vector F. To illustrate, consider a network with two OD
pairs and having three paths pi, p2' p3 joining the first pair and four
paths p 4 p', p6 and p7 joining the second pair, where the pi are of
84
dimension R. These give rise to 12 (3 times 4) aggregate extreme
flows. Now let F = .5p1 + .5P2 + Op3 + .3p4 + .2p5 + .5p6 + Op7 be a
flow pattern expressed in terms of five path flows. Since F = .3(pl +
p4) + .5(p2 + p6), the same flow pattern is expressible in terms of
just three aggregate extreme flow vectors (pl + 4 ), (p + p5), and
(P2 + P6). Thus working in the aggregate extreme flow space can save
computer storage, a paramount concern on large networks.
Another difference in the algorithms is that Bertsekas and
Gafni only do one projection iteration (in the path flow space) before
generation of new paths (via a subproblem which is equivalent to Step 1
of our algorithm). Our Step 2 does as many projection iterations as
necessary to achieve an et solution. This tradeoff takes advantage of
the fact that on large networks the generation of shortest paths in
Step 1 is far more time consuming than the projection iterations.
The computational experiments summarized in the next section
show that even on small networks there are computational savings in
both storage and computer time.
4.2.3 Computational Testing
In this section, computational testing of the algorithm on
four small problems is summarized. In the proof of convergence, we
have assumed that et + 0 and that the stopping rule is G(xt) = 0. For
implementation purposes it is more practical to require an esolution.
Accordingly, we chose et = 10 + 5 104/t and 6 = 104. Further
more, the algorithm is terminated when the value of the relative gap,
t t t t t 4
defined to be C(x)(x yt)/C(x )y is less than 104. This is a
more appealing choice than G(xt) and is further justified by the re
sults for problem ND1 (below). In all testing we solve the master
85
variational inequality problem in Step 2 by the projection method. The
matrix which defines the norm for the projection mapping does not
change from iteration to iteration, and it is taken to be a diagonal
matrix whose diagonal elements have the same positive value.
All programs are written in FORTRAN and were executed on the
IBM 3033N computer at the University of Florida. In order to make an
approximate comparison with the times reported by Nguyen and Dupuis
[1981], we must relate their times on a Cyber 173 machine to the IBM
3033N. The ratio of the machine cycle times for the IBM 3033N and the
Cyber 173 is 58/50 = 1.16. However, when we solved their symmetric
traffic equilibrium by the program TRAFFIC written by Nguyen and James
[1975] on the IBM 3033N, we found that the average ratio of the time on
the IBM 3033N and the Cyber 173 is approximately 3. We cannot make any
definite conclusion from these facts except that IBM 3033N is clearly a
slower machine.
Two algorithms are selected for the comparison, the algo
rithm by Nguyen and Dupuis [1981] and a variation of the algorithm pre
sented in Bertsekas and Gafni [1982J. Other techniques, which have as
their subproblem a symmetric traffic assignment problem, require a stand
ard nonlinear programming algorithm, e.g., the FrankWolfe algorithm,
for its solution. Such techniques would construct many shortest path
trees per iteration and therefore are likely to be too time consuming
for a large network. The reader interested in these techniques is re
ferred to the paper by Fisk and Nguyen [1982].
The algorithm by Nguyen and Dupuis [1981] is an interesting
modification of the algorithm for nperson games by Zuhovicki et al.
[1969]. It employs the cuttingplane algorithm to solve the optimi
zation version of the variational inequality problem, VIP[C,S], i.e.,
86
max min C(x)(x u)
ucS xeS
A relaxed version of this maxmin problem can be written as a linear
programming problem
max x
subject to
C(x (xi U) 2 x i = 1, 2, ..., n
U C S,
i
where x i = 1,..., n, are points belonging to S. During each iter
ation, such a relaxed problem is solved by the simplex method, and the
resulting solution is used to generate another constraint (a cut),
n+l n+l
C(xn+ )(xn u) 2 x0,
to be included in the next relaxed linear programming problem. In
t
their computational testing x0, the value of the objective function at
iteration t, monitors the progress of the algorithm.
For comparison purposes we have also implemented a variation
of the Bertsekas and Gafni method which we refer to as the commodity
projection method. (As previously described, their method works in the
path flow space and does not involve column dropping.) In the commod
ity projection method, a flow vector is retained for each "commodity"
which is defined to be the aggregate flow from an origin to all desti
nations. This strategy only requires a projection for each origin
at each iteration, rather than for each origindestination (OD) pair.
We only report here the results for the "allatonce" version proposed
by Bertsekas and Gafni because it performed better than the "oneata
time" version (see the discussion in Bertsekas and Gafni [1982]) on our
87
test problems. In implementing the commodity projection method the
choice of the projection matrix is as described before.
The first 3 problems, ND1, ND2, and ND3, are taken from
Nguyen and Dupuis [1981]; all of them have the same network topology
and trip matrix as shown in Figure 41. This network admits 8 paths
from node 1 to node 3, 6 paths from 1 to 4, 4 paths from 2 to 3 and 6
paths from 2 to 4. Thus, there are a total of 8.6.4.6 = 768 possible
extreme flow patterns. The cost vectors for ND1, ND2 and ND3 corres
pond to the symmetric, diagonally dominant and nondiagonally dominant
cost vectors in Nguyen and Dupuis [1981].
In Problem ND1 the Jacobian C(x) is diagonal; therefore an
equivalent optimization problem exists. Table 43 summarizes the com
putation results for the simplicial algorithm as compared with the re
sults reported by Nguyen and Dupuis. Judging from the objective val
ues, the two methods produce equivalent solutions, with thesimplicial
algorithm time offering an improvement even though the IBM 3033N is the
slower machine. Of course this comparison is tentative because the
computers differ and because both algorithms can be improved by "fine
tuning." In particular, the Nguyen and Dupuis method could be modified
to solve the linear program approximately at each iteration and this
might improve its overall performance.
Note, also, that the relative gap provides a more accurate
measure of the solution than the gap itself.
Tables 44 and 45 summarize the results for Problem ND2 and
ND3, respectively. Note that out of 768 extreme flow patterns the
simplicial algorithm generates no more than 6 extreme flow patterns and
retains at most 4 patterns during any iteration. For problems ND2, the
Table 43: Computational Results for Problem ND1
Nguyen and Oupuis The Simplicial Algorithm'
t F(xt) (CPUmele.) t F(xt) Gap Rel. Gap p p (CTimec.)
1 0 1 CPUsec') 11 2 (CPU sec.)
100,400.00
25,285.43
7,189.15
2,451.38
2,284.28
669.05
155.94
27.48
3.90
.90
.29
94,253.37
86,847.30
86,294.06
86,082.69
85,963.00
85,210.60
85,061.23
85,035.67
85,030.29
85,028.53
85,028.14
0.016
0.041
0.071
0.093
0.120
0.334
0.495
0.640
0.750
0.862
0.975
94,252.81
86,560.56
85,894.68
85,513.25
85,214.62
85,046.81
85,027.75
85,027.62
85,027.62
.1004 10"
.2528 105
.4844 104
.5183 104
.1835 104
.1849 104
.6565 103
.6744 102
.6750 101
.1136 101
.2823 100
.4951 101
.5343 101
.1855 10'1
.1873 101
.6558 102
.6711 103
.6713 104
.0910
.1265
.1722
.2164
.2571
.3313
.4263
.4887
.5178
Notes: 1: Time
2: Pl
3: p2
on Cyber 173
I of active extreme aggregate
# of projections
flow patterns
4: Time on IBM 3033
5: 150 is the maximum number of projections allowed per iteration
Table 44: Computational Results for Problem ND2
Nguyen and Oupuls The Simplicial Algorithm The Commodity Projection Method
t o Tie t Re. ap. P P Time2 t Rel. Gap. p5 Time2
0 Time t Rel. Gap. p1 p2 Time2 t IR. Gap. p3
.27 in1 .02
2
3
3
4
4
.0882
.1209
.1534
.2038
.2329
+
Notes: 1: Time on Cyber 173
2: Time on IBM 3033
3: pl a I of active extreme aggregate flow patterns
4: p2 
5: p3 I
.2677 101
.2035 101
.5091 10'2
.1449 10.2
.4229 103
.1241 103
.9681 104
.9681 .7111
of projections
of active extreme commodity
flows for commodities 1 and 2
4,2
4,2
4,2
4,2
4,2
4,2
.2677 101
.3446 100
.1169 100
.1852 101
.6524 104
176,492.00
22,448.71
12,156.66
3,666.34
2,685.25
1,511.72
972.88
715.99
195.13
133.25
87.26
71.98
30.82
25.58
0.03
0.01
.026
.046
.070
.101
.127
.146
.168
.190
.222
.255
.278
.300
.322
.342
.364
.387
.0827
.1847
.2846
.4355
.5599
.6862
.7111
Table 45: Computational Results for Problem ND3
Nguyen and Dupuis The Simplicial Algorithm The Commodity Projection Method
t x Time1 t Rel. Gap. p p2 Time2 t Rel. Gap. p Time2
0 P3 t 3
113,420.00
14,909.32
1,731.35
773.63
373.64
285.58
178.94
113.39
90.73
26.94
8.23
0.76
0.70
0.40
0.33
0.03
0.01
0.023
0.046
0.068
0.091
0.111
0.138
0.159
0.198
0.238
0.258
0.282
0.313
0.336
0.359
0.380
0.401
0.430
.1941
.4185
.8564
.6111
.1637
.6401
89
0.0856
0.1181
0.1567
0.1963
0.2318
0.2604
.1941 11n
.6743 101
.6768 102
.2409 102
.9164 103
.3479 103
.1326 10"3
.9026 104
0.0807
0.1801
0.3052
0.4306
0.5571
0.6827
0.7820
0.8575
on Cyber 173
on IBM 3033
I of active extreme aggregate flow patterns
: p I # of projections
: P3 of active extreme commodity
flows for commodities 1 and 2
Notes: 1: Time
2: Time
3: p, "
91
simplicial algorithm built 6 minimum path trees and performed 109 pro
jections, whereas the commodity projection method built 26 minimum path
trees and performed 26 projections. Although the simplicial algorithm
performed many more projections, the CPU time indicates that the time
required to build the additional 20 minimum path trees more than off
sets the time to do the extra projections. A similar conclusion can be
drawn from Table 45. The computation times of the simplicial algo
rithm are the best on both ND2 and ND3, but, for the reasons already
given, comparison with the Nguyen and Dupuis method must be regarded as
tentative.
The network topology and the trip matrix for Problem 4 are
shown in Figure 47. The network is from Fisk and Nguyen [1982], but
we constructed the cost function shown in Table 46. This network ad
mits 4 paths between each OD pair with nonzero demand; in other words,
there are (4)12 = 16,777,216 possible extreme flow patterns. Table 47
summarizes the computational results of the simplicial algorithm versus
the commodity projection method. The conclusions from Problem ND2 and
ND3 also hold for this case. Considering the very large number of ex
treme flow patterns it is encouraging that the simplicial algorithm
generated only 18 extreme flow patterns.
We feel the above results demonstrate the potential of the
simplicial algorithm with regard to computer time, and they show fur
ther that the method obtains good solutions while generating only a
small number of extreme flow patterns.
92
Network Topology
5 6 7 8
1 0 35 30 20
2 40 0 35 25
3 30 20 0 35
4 35 25 40 0
Trip Matrix
Figure 47: Data for Problem 4
Table 46: Cost Vector for Problem 4
Arc 0 0 Ci(x) Arc 0 0 C (x)
1 1 9 f(x1 + 14) 19 14 11 f(1/2 x16 + 19 + x25 + 1/4 x26)
2 1 12 f(x2 + 1/2 x11 + 1/4 x12 + 1/4 x14) 20 14 18 f(x16 + x20)
3 2 10 f(1/4 x10 + x3 + x21 + 1/4 x22) 21 15 10 f(x3 x21)
4 2 16 f(xo0+ x4) 22 15 16 f(1/4 x10 + 1/4 x3 + x2)
5 3 19 f(x34 + 1/4 x33 + 1/4 x32 + x5) 23 16 14 f(x29 + 1/4 x30 + 1/4 x18 + x17 + x23)
6 3 22 f(x32 + x6) 24 16 20 f(x18 + x24)
7 4 17 f(x7 + x36) 25 17 11 f(x19 + x25)
8 4 21 f(x8 + 1/2 x27 + 1/4 x28 + 1/4 x36) 26 17 13 f(x15 + 1/4 x16 + 1/4 x19 + x20 + x26)
9 9 6 f(xg + x22) 27 18 8 f(1/4 x8 + x27 + x35 + 1/4 x36)
10 9 16 f(x10 + 1/2 x3 + 1/4 x4 + 1/4 x22) 28 18 21 f(xg + x28)
11 10 5 f(1/4 x2 + x11 + x13 + 1/4 x14) 29 19 14 f(x29 + x23)
12 10 12 f(x2 + x12) 30 19 15 f(x30 + 1/4 x18 + 1/2 x17 + 1/4 x23)
13 11 5 f(x11 + x13) 31 20 7 f(x33 + x31)
14 11 9 f(xI + 1/4 x2 + 1/4 xll + x14) 32 20 22 f(1/4 x33 + x32 + x6 + 1/4 x5)
15 12 13 f(x15 + x26) 33 21 7 f(x33 + 1/2 x32 + 1/4 x31 +1/4 xS)
16 12 18 f(x16 + 1/4 x19 + 1/2 x20 + 1/4 x26) 34 21 19 f(x34 + x5)
17 13 15 f(x30 + x17) 35 22 8 f(x27 + x35)
18 13 20 f(1/2 x30 + x18 + x24 + 1/4 x23) 36 22 17 f(x7 + 1/4 x8 + 1/4 x27 + x36)
f(Z) 3.0 (1 + 0.15(Z/35)4)
Table 47: Computational Results for Problem 4
The Simpllcial Algorithm The Commodity Projection Method
t Rel. Gap. p2 p2 Time t Rel. Gap. p3 Time1'
2
3
4
5
6
5
5
6
7
7
8
9
10
8
9
10
11
12
12
23
29
85
82
78
161
208
3005
174
3005
3005
3005
193
3005
169
249
186
180
0.1092
0.1543
0.2440
0.3409
0.4434
0.6078
0.8118
1.1308
1.3477
1.7025
2.0920
2.5201
2.8264
3.2162
3.4686
3.8587
4.1790
4.5118
4.5390
Notes: 1: Time on IBM 3033
2: pl I of active extreme aggregate flow patterns
3: P2 I of projections
.3624
.2717
.1957
.1405
.1014
.7753
.6118
.5441
.4896
.4351
.3908
.3539
.3227
.3008
.2829
.2667
.2522
.2390
.2269
.2159
.2057
2,2,2,2
2,3,2,2
2,3,2,2
2,3,3,3
2,3,3,3
2,4,3,3
2,4,3,3
2,5,3,3
2,5,3,3
2,5,3,3
2,5,3,3
2,5,3,3
2,5,3,3
2,6,3,3
2,6,3,3
2,6,3,3
2,6,3,3
2,5,3,3
2,6,3,3
2,6,3,3
2.6.3.3
.1007
.3587
.6445
.9286
1.2188
1.5054
1.7952
2.0835
2.3732
2.6625
2.9534
3.2440
3.5347
3.8285
4.1208
4.4118
4.4018
4.9910
5.2798
5.5692
5.8599
4: p3 I of extreme commodity flows
5: 300 max. allowable number
of projections
.3624
.1018
.4250
.2450
.1952
.9912
.1489
.1082
.1222
.2853
.1133
.4488
.4447
.1071
.1116
.4216
.1822
.1498
.6774
3317
CHAPTER FIVE
SUMMARY AND FUTURE RESEARCH
5.1 Summary of Results
This study is concerned with decomposition techniques for the
variational inequality formulation of the traffic assignment problem.
This formulation includes the standard (optimization) formulation which
has been investigated previously by many authors. In this disserta
tion, it is shown that two decomposition techniques, transfer and sim
plicial decompositions, originally developed for the optimization
formulation, can be extended to the variational inequality formulation.
Chapter 3 discusses the extension of transfer decomposition
to variational inequality problems. First, the modified version of
transfer decomposition for the variational inequality problem is pre
sented, and it is shown that transfer decomposition provides a tech
nique to separate the symmetric component from the rest of the problem.
Then, several properties for the mapping which defines the master vari
ational inequality problem are stated and proved. The main property is
that the mapping resulting from applying transfer decomposition to the
traffic assignment problem is always pointtopoint. Using this prop
erty, an existing algorithm can be applied to the master variational
inequality. In the last section of this chapter, it is shown that
transfer decomposition also provides a rigorous basis for constructing
aggregate functions used in network simplification.
95
