To be safe, one could approach problem (36) via some tree- search scheme. This would surely yield the exact optimal solu- tion, but the computing time required might become prohibitive.
On the other hand, one could try to develop some heuristic approaches exploiting submodularity as far as possible, hoping to find at least a good local maximum. This would usually
require negligible computing time, but unfortunately submodularity alone does not provide any sufficient condition for a global
maximum.
However, the successful experience with simple heuristics applied to maximizing submodular functions (although of a less general nature than the one considered here, see Wolsey, 1980) suggests that the actual performance of such heuristics could be
very good. In an exploratory stage of research, it is therefore worth checking the results obtained with heuristic algorithms against the result obtained with an exact method. If numerical experience shows that heuristics work almost or just as good as the exact algorithms, then good reasons for further theoretical investigation are provided.
A tree search can be easily organized by exploiting proposi- tions 7-11. Suppose, for instance, the search goes downward,
i.e., building smaller subsets. Then the search within a given subset S can be stopped when:
a. S is a local maximun;
b. S is a nondecreasing point (see Proposition 10);
c. an upper bound to G(T), T
-
C Sf as computed by (41), is less than the highest value of the objective function found so far.For condition c, an upper bound can be computed from (41 ) as follows
for all T - C S
where
The procedure outlined above is based on the submodularity property only. Another method should be mentioned, which has
been developed and used in Erlenkotter and Leonardi (forthcoming).
This method, referred to as INTLOC, does not use any submodularity at all, but rather works with the continuous relaxation of (36) to get approximate integer solutions and bounds in the tree search.
More precisely, let the function 9(X) = E P i $ i ( ~ + l o g ~ )
-
E x . ai j I j
be defined over all real nonnegative vectors X = x , where the
4 . (V) 1 are the functions defined by equation (1 9)
.
Then, ifit follows from (30), (31), and (32) that
Therefore the solution to the mathematical program max .(R(X) : O
-
< x j - < 1 ,Vj)X
provides an upper bound to the optimal solution of
If problem (48) is solved by a simple Frank-Wolfe method, feasible integer solutions are also generated as a by-product at each iteration. The best of such solutions can therefore be taken as an approximation to the solution of the problem
max CQ(X) : x E {0,1)
,
~ j )X j
which is of course equivalent to (36).
Such a procedure can be easily embedded into a branch-and- bound search scheme. In its present version INTLOC is not
designed to work with a general function of type (46). It
assumes a special structure, the same one described in Section 6 of this paper.
The heuristic procedures will be kept as simple as possible.
The first proposed heuristic tries to find a weak local maximum, while monotonically increasing the value of the objective function.
This can be done as follows. Assume a given iteration starts with
a t r i a l s o l u t i o n s e t S which d o e s n o t meet c o n d i t i o n s ( 4 2 ) . T h i s means t h a t t h e c u r r e n t v a l u e o f t h e o b j e c t i v e f u n c t i o n , G ( S ) , c a n b e i n c r e a s e d by a d d i n g o r d r o p p i n g o n e e l e m e n t . L e t t h e c h a n g e g i v i n g t h e h i g h e s t i n c r e a s e i n t h e o b j e c t i v e f u n c t i o n b e i n t r o - d u c e d and r e p e a t t h e s t e p . The p r o c e d u r e s t o p s when n o e l e m e n t c a n b e c o n v e n i e n t l y added o r d e l e t e d , t h a t i s , when a weak l o c a l maximum i s d e t e c t e d . T h i s p r o c e d u r e w i l l b e s t a t e d f o r m a l l y :
HEURISTIC 1 . ( A s c e n t t o w a r d s a weak l o c a l maximum) 1 . g u e s s a s t a r t i n g S
2. f i n d p i ( S ) = max p
.
( S ) andjgs
Ip k ( S
-
{ k l ) = min p . (S-
{ j l )j E s I
3 . i f p i ( S )
-
< 0 a n d p k ( S-
{ k ) )2
0 s t o pi f p i ( S ) >
-
p ( s - { k } ) r e p l a c e S by s U { i } ki f pi ( S )
-
< - p k ( S-
{ k ] ) r e p l a c e S by S-
{ k)The s t a r t i n g g u e s s i s q u i t e a r b i t r a r y , and d i f f e r e n t s t a r t s may l e a d t o d i f f e r e n t l o c a l maxima. When no b e t t e r s t a r t i s a v a i l - a b l e , a r e a s o n a b l e o n e i s
The s e c o n d h e u r i s t i c t r i e s t o f i n d a s t r o n g l o c a l maximum, w h i l e m o n o t o n i c a l l y i n c r e a s i n g t h e v a l u e o f t h e o b j e c t i v e f u n c t i o n . I t works a s f o l l o w s . Suppose a g i v e n i t e r a t i o n s t a r t s w i t h a weak l o c a l maximum. I f a t l e a s t o n e p a i r e d i n t e r c h a n g e i m p r o v e s t h e
v a l u e o f t h e o b j e c t i v e f u n c t i o n , a new and b e t t e r weak l o c a l maximum i s produced by r e s t a r t i n g H e u r i s t i c 1 w i t h t h e i n t e r - changed s o l u t i o n , and t h e i t e r a t i o n i s r e p e a t e d . The p r o c e d u r e s t o p s when no p a i r e d i n t e r c h a n g e can improve t h e c u r r e n t weak l o c a l maximum, which w i l l t h e r e f o r e b e a s t r o n g l o c a l maximum t o o .
T h i s p r o c e d u r e w i l l be s t a t e d f o r m a l l y :
HEURISTIC 2 . ( A s c e n t t o w a r d s a s t r o n g l o c a l maximum)
1 . u s e H e u r i s t i c 1 w i t h any s t a r t t o g e n e r a t e a weak l o c a l maximum S
2. s e t So = S and T = S 3. i f T = % s t o p
4 . choose some j E T and r e p l a c e T by T
-
{ j )5. u s e H e u r i s t i c 1 w i t h s t a r t S o - { j ) t o g e n e r a t e a weak l o c a l maximum S
J u s t a s f o r H e u r s i t i c 1 , no u n i q u e n e s s o f t h e f i n a l s o l u t i o n i s a s s u r e d i n p r i n c i p l e . However, a l a t e r s e c t i o n on n u m e r i c a l e x p e r i m e n t s w i l l show t h a t t h e performance o f H e u r i s t i c 2 i s s u r p r i s i n g l y b e t t e r t h a n m i g h t be f o r e s e e n from t h e o r y .
A SPECIAL CASE
The r a n d o m - u t i l i t y model c o n s i d e r e d s o f a r i s q u i t e
g e n e r a l , s i n c e no r e s t r i c t i v e assumption h a s been i n t r o d u c e d f o r t h e d i s t r i b u t i o n f u n c t i o n s F ( Y )
.
A simplifying assumption often found in the literature (Domencich and McFadden, 1975; Daly, 1978; and Van Lierop and Nijkamp, 1979) is that Y is a sequence of independent identically distributed random variables with a common distribution function
where a and B are called the shape and the location parameters, respectively. This distribution is known as the Gumbel distri- bution, and it plays an important role in extreme order statistics
(Gumbel, 1958; and Galambos, 1978).
The mean of (5 1 ) is known to be IJ = xdD (x) = - B
+
ya a
where y = 0.5772157... is Euler's constant.
Because of independency, F(Y) takes the form B -aYj
F(YI = IID(yj) = exp [-e i e
j j
I
and the extreme value distribution is, according to (18)
F (x
-
V) = exp [-e' h(v) e-ax] = exp 1-e-
[ax-
B-
109 h ( ~ )I I
(53) whereComparison of (53) with (51) shows that F(x-V) is still a Gurnbel distribution, with shape parameter a and location para- meter B + log h(V). Therefore, according to (19) and (52), the expected utility for a customer in a given origin is
0 (V) = xdF (x
-
V) =-
a 1 log h(V)+ -
a B+ v
- aSince, because of Propositions 1 and 2, a shift in the origin of the utility scale does not affect customer choice behavior, addi- tive constants can be dropped from (55) and the expected utility can be redefined as
The choice probabilities can be found using (23) and (54). They are
This is the well known multinominal logit model, extensively used in transport demand analysis (Domencich and McFadden, 1975;
Williams, 1977).
Let all these results be introduced in problem (36). From (30) one gets
Substitution in (3 1 ) yields
1 avi j
g(L) = ,E Pi log E e i j EL
and substitution of (58) in (32) gives the objective function G(L) =
c , ~
1 pi log E eavij-
L ai j EL j a j
Since no generality is lost if the fixed charges are rescaled by a and G(L) is multiplied by a, the objective function can be
redefined as
where
Under the above assumptionsl problem (36) takes the form max C Pi log L fij
-
C aL i j EL jEL 1
Problem (60) may be given an alternative formulation, in the spirit of ~roposition 5. Since this formulation is closely rela- ted to the one extensively used by Coelho (1977, 1979, 1980a, and 1980b) the following will be called
PROPOSITION 12. (Coelho representation)
P r o o f
the proof parallels the one of Proposition 5.
In order to implement the algorithms of Section 5, compu- tational formulas for the incremental values p.(S) and p.(S
-
(j))3 3
a r e needed. From ( 3 7 ) and ( 5 9 ) t h e s e a r e
7. SOME NUMERICAL RESULTS
The e x a c t and h e u r i s t i c a l g o r i t h m s o f S e c t i o n 5 h a v e b e e n a p p l i e d t o a t e s t problem. The problem d a t a r e f e r t o t h e l o c a t i o n o f h i g h s c h o o l s i n T u r i n , I t a l y . A d e t a i l e d d e s c r i p t i o n o f t h e d a t a and t h e g e o g r a p h i c a l s e t t i n g c a n b e f o u n d i n L e o n a r d i and B e r t u g l i a (1981) and i n E r m o l i e v , L e o n a r d i , a n d V i r a ( 1 9 8 1 ) , where t h e y h a v e b e e n u s e d t o t e s t somewhat d i f f e r e n t o p t i m a l
l o c a t i o n a l g o r i t h m s (namely, a problem w i t h c o n s t r a i n t s on
f a c i l i t y s i z e and a s t o c h a s t i c programming a p p r o a c h ) . The i n p u t d a t a u s e d i n t h e t e s t s a r e r e p o r t e d i n t h e Appendix. From t h e p o i n t o f v i e w o f n u m e r i c a l t e s t i n g , s a l i e n t f e a t u r e s o f t h e problem a r e :
a . T u r i n i s d i v i d e d i n t o 23 d i s t r i c t s , e a c h d i s t r i c t b e i n g b o t h a p l a c e o f r e s i d e n c e o f h i g h s c h o o l demand and a p o s s i b l e l o c a t i o n f o r a h i g h s c h o o l f a c i l i t y . No l i m i t a t i o n i s p l a c e d on t h e c h o i c e o f d i s t r i c t s o f
d e s t i n a t i o n , s o t h a t i n p r i n c i p l e a c u s t o m e r l i v i n g i n a g i v e n d i s t r i c t c a n u s e a f a c i l i t y i n any o t h e r d i s t r i c t . b. The u t i l i t y h a s b e e n s i m p l y s e t e q u a l t o t r a v e l t i m e
changed i n s i g n , i . e . :
where
C ij is the travel time from district i to district j, in minutes
Travel times are measured on the public transport net- work. The within-district travel time has been given a
standard value of five minutes, in accordance with empirical findings.
The parameter a has been given a value
According to more recent origin-destination surveys on home-to-school trips, this parameter should be set equal to values around 0.15.
However, the value (65) has been kept, in order to make results comparable with previous studies (Erlenkotter and Leonardi, forthcoming).
c. The quantities Pi appearing in (60) are the number of high school students living in each district, as of 1977
(Provincia di Torino, 1978). The values of these quan- tities range from 500 to 2,500, approximately.
d. The fixed charges a have been set equal for all j
districts
and the algorithms have been tested for values ranging from 500 to 5,000. his range has been used for testing purposes only [the difficulty of (60) usually increases with a], and there is no claim of realism in it.
The r e s u l t s o f t h e n u m e r i c a l t e s t i n g a r e summarized i n
Table 1. Comparison of the performance of exact and heuristic algorithms for the Turin high school. test problem.
O b j e c t i v e f u n c t i o n (changed i n s i g n ) CPU time (seconds) a
F i r s t
Fixed Tree s t a g e of Heuris- Heuris- Tree C Heuris- Heuris-
charge s e a r c h INTLOC t i c 1 t i c 2 s e a r c h t i c 1 tis 2
500 25899 25899 25899 2 5899 1.2 4.1 15.7
a On t h e IIASA VAX computer.
bBoth w i t h t h e s e a r c h based on submodularity bounds and w i t h INTILX: [ t h e Frank-Wolfe based branch-and-bound a l g o r i t h m used i n E r l e n k o t t e r and Leonardi ( f o r t h c o m i n g ) ] ; t h e r e s u l t s a r e t h e same, and correspond t o t h e e x a c t s o l u t i o n .
C Recorded CPU times r e f e r t o t h e s e a r c h based on submodularity bounds, which p r o v i d e d an answer i n r e a s o n a b l e t i m e f o r f i x e d c h a r g e s up t o 3000. For h i g h e r f i x e d charge v a l u e s t h e computing t i m e was t o o long t o w a i t f o r an answer. Although e x a c t t i m e s f o r INTLOC were n o t r e c o r d e d , i t had t o be run one n i g h t t o p r o v i d e t h e r e s u l t s f o r h i g h f i x e d charge values.
previously produced with INTLOC (which was just as slow, although its users were less impatient) was kept.
The CPU time of tree search methods is very low for small fixed charge values, but it starts a fast increase after a fixed charge of 2,500 and becomes infeasible for fixed charge values higher than 3,000.
The CPU times for Heuristic 1 are all around 2 - 4 seconds, no matter what the value of the fixed charge. (In Table 1, they seem to decrease with the fixed charge, but this cannot be generalized, since it depends on the starting solution used.) Heuristic 1 seems therefore a bit slower than the tree search for small fixed charge values, but this is more than counter- balanced by its performance for high fixed charge values.
The CPU times for ~euristic 2 vary roughly between 10 and 20 seconds, again independently of the fixed charge value. They are not negligible, but they are still much lower than the ones for the tree-search algorithms in the hard cases. This must be
coupled with the fact that apparently ~euristic 2 n e v e r fails to reach the optimum,
Although the numerical tests discussed above cannot be claimed to be exhaustive, they seem to be enough to state that
a. Heuristics 1 and 2 provide a uniformly superior start than any other method for a further tree-search
refinement.
b. The high performance of Heuristic 2 deserves further theoretical investigation, along the lines of looking for possible sufficiency conditions.
c. For real sensitivity analysis problems, Heuristic 2 can be safely recommended; for a preliminary rough analysis, Heurisitc 1 can be enough,
Of course the validity of the above statements is provisionally limited to objective functions of the form (60), although it might be argued that the performance would be just as good with many other submodular set functions.
Statement a is perhaps the most intriguing one. Indeed the tree search starting with the results of Heuristic 2 has been attempted, but this has led to no significant improvement in its performance. In other words, recognizing the starting solution as an optimal one seems to take just as long as building an
optimal solution from a bad start. Improving the tree search and tightening its bounds is therefore another subject for further investigation.
Statement b is related to the possible development of an effective duality theory for problems of type (60), or more generally of type (8). The effectiveness of dual relationships is the main reason for the successful algorithms recently devel- oped to solve problems of type (1)-(3) (Bilde and Krarup, 1977;
Erlenkotter, 1978), and this encourages the search for similar results for the more difficult problems (8) and (60).
Statement c reminds us that the mathematical interest in finding an exact solution is not necessarily realistic. Assuming an impatient decision maker would have used only'Heuristic 1, according to Table 1 his maximum loss (in terms of relative difference between obtained and optimal objective value) would have been (for the fixed charge value of 2,500)
that is about 0.3%. How many input data have measuring errors less than this figure?