The Use of Random-Utility Theory in Building Location-Allocation Models

NOT FOR QUOTATION WITHOUT PERMISSION OF THE AUTHOR

THE USE OF RANDOM-UTILITY THEORY I N BUILDING LOCATION-ALLOCATION MODELS

*

G i o r g i o ~ e o n a r d i March 1981

WP-81-28

To b e p r e s e n t e d a t t h e 9 t h T r i e n n i a l C o n f e r e n c e on O p e r a t i o n a l R e s e a r c h , I n t e r n a t i o n a l F e d e r a t i o n o f Opera- t i o n a l R e s e a r c h S t u d i e s (IFORS ' 8 1 ) , Hamburg, FRG, J u l y 20-24, 1981

*Working P a p e r s u b m i t t e d f o r p u b l i c a t i o n i n t h e f o r t h c o m i n g book S p a t i a l A n a l y s i s of Pub l i c F a c i l i t i e s , e d i t e d by J . T h i s s e and H. Z o l l e r . Working Papers a r e i n t e r i m r e p o r t s on work o f t h e I n t e r n a t i o n a l I n s t i t u t e f o r A p p l i e d S y s t e m s A n a l y s i s and have r e c e i v e d o n l y l i m i t e d r e v i e w . V i e w s o r o p i n i o n s e x p r e s s e d h e r e i n d o n o t n e c e s s a r i l y r e p r e - s e n t t h o s e o f t h e I n s t i t u t e o r o f i t s N a t i o n a l Member O r g a n i z a t i o n s .

INTERNATIONAL INSTITUTE FOR APPLIED SYSTEMS ANALYSIS A-2361 Laxenburg, A u s t r i a

THE AUTHOR

D r . G i o r g i o L e o n a r d i h a s b e e n a t IIASA s i n c e O c t o b e r 1979 on l e a v e from t h e P o l y t e c h n i c a l I n s t i t u t e o f T u r i n , I t a l y . H e h a s p r e v i o u s l y b e e n a r e s e a r c h e r a t t h e I t a l i a n A s s o c i a t i o n f o r Housing R e s e a r c h (AIRE) a n d a c o n s u l t a n t on u r b a n a n d r e g i o n a l p l a n n i n g f o r t h e R e g i o n a l A u t h o r i t y o f P i e m o n t e .

FOREWORD

The public provision of urban facilities and services often takes the form of a few central supply points serving a large number of spatially dispersed demand points: for example, hospitals, schools, libraries, and emergency services such as fire and police. A fundamental characteristic of such systems is the spatial separation between suppliers and consumers. No market signals exist to identify efficient and inefficient geo- graphical arrangements, thus the location problem is one that arises in both East and West, in planned and in market economies.

This problem is being studied at IIASA by the Public Facility Location Task, which started in 1979. The expected results of this Task are a comprehensive state-of-the-art survey of current theories and applications, an established network of international contacts among scholars and institutions in different countries, a framework for comparison, unification, and generalization of exist- ing approaches, as well as the formulation of new problems and

approaches in the field of optimal location theory.

This paper is both a unifying effort and a contribution to the formulation of such new problems and approaches. It explores the relationships between a recent area of geographic research, random-utility theory, and a recent area of applied mathematic research, the optimization of submodular functions. The frdit- fulness of this marriage is shown by some numerical results, which seem to suggest that the approach can yield new, powerful tools for location problems.

Related publications in the Public Facility Location Task are listed at the end of this report.

Andrei Rogers Chairman

Human Settlements and Services Area

ABSTPACT

The most important part of a location-allocation model is the allocation rule, that is, the way clients are assigned to facilities. In the well-known models of the "plant-location"

family, the embedded allocation rule is the assignment of the least-travel-cost facility,

This allocation rule depends on the assumption that the cost, or more generally utility, associated with each possible facility choice is deterministically known. The simplest way to generalize a plant-location model is to add a random term to travel costs, with a known probability distribution. Such

randomness may be shown to arise in many real-life situations, and the resulting choice models constitute the subject of

random-utility theory,

This paper introduces the use of the random-utility modeling philosophy in location-allocation problems, Some relevant prop- erties of the resulting family of models are derived, Among them, of special importance is the submodularity property, which relates the random-utility-based location models to a recent area of research in combinatorial optimization. Submodularity is

exploited to develop simple heuristic algorithms, and the effec- tiveness of the approach is supported with some.numerica1 results.

CONTENTS

1. INTRODUCTION, 1

2. A GENERALIZATION OF THE UNCAPACITATED FACILITY LOCATION PROBLEM, 2

3. BASIC PROPERTIES OF THE ADDITIVE RANDOM-UTILITY MODEL, 6 4. BASIC PROPERTIES OF SUBMODULAR SET FUNCTIONS, 15

5. SOME ALGORITHMS, 21 6. A SPECIAL CASE, 25

7. SOME NUMERICAL RESULTS, 29

8. CONCLUDING COMMENTS AND ISSUES FOR FUTURE RESEARCH, 34 APPENDIX, 36

REFERENCES, 37

THE USE OF RANDOM-UTILITY THEORY IN BUILDING LOCATION-ALLOCATION MODELS

1. INTRODUCTION

Some recent theoretical and computational innovations have revived interest in location problems in the last few years.

On the theoretical side, the need to introduce more realis- tic and economically sound measures of customer benefits has been felt by many scholars. An outstanding contribution in this field has been given by the use of random-utility theory. This approach, mainly developed for transport demand analysis

(Domencich and McFadden, 1975; ~illiams, 1977), has been rapidly extended to the spatial allocation of economic activities

(~oelho,' 1977, 1979; Coelho and Williams, 1978; ~ a c ~ i l l and ~ilson, 1979). Its use in the facility location context first appeared

in Coelho and Wilson (1 976) and independently in Leonardi (1975)

,

then, shortly after in many other contributions (Williams and Senior, 1977; Leonardi, 1978, 1980a, 1980b; Beaumont, 1979, 1980;

Coelho, 1980a)

.

On the computational side, the major event has been the

appearance of the dual ascent method, which amazingly outperforms all previously used methods to solve the classic uncapacitated plant-location problem. The method, first suggested by Bilde and Krarup (1977), has been developed by Erlenkotter (1978) and

Van Roy and Erlenkotter (1980). More interesting than the method itself have been the theoretical investigations it stimulated, Among them, the most fruitful one has been the parallel develop- ment of a new framework to analyze cornbinatorial optimization

problems with submodular objective functions. The method has been developed by Cornuejols, Fisher, and Nemhauser (1977); Nemhauser, Wolsey, and Fisher (1978); Fisher, Nemhauser, and Wolsey (1978);

and its relationships with the"dua1-+scent method have been analyzed in Wolsey (1980).

This paper tries a first step in putting the two theories together, Submodularity (the key property leading to the success of the above methods) is shown to hold for all objective functions based on the additive random-utility model. Some simple heuris- tics based on this property are proposed, and some (surprisingly good) numerical results are shown for a typical random-utility uncapacitated location problem.

As a conclusion, it is argued that further investigation could yield a substantial improvement in the state-of-the-art of location mo6eling for the near future, and provide superior

algorithms for a much wider class of problems than the ones usually considered in the operations research literature.

2. A GENERALIZATION OF THE UNCAPACITATED FACILITY LOCATION PROBLEM

The classic uncapacitated facility location problem

(Efroyrnson and Ray, 1966; Spielberg, 1969) can be formulated in the following way

min C C xij Cij

+

C a L,X j f ~ i j EL j

where

i labels customer locations j labels facility locations

Pi is the total number of customers in location i x is the number of customers in location i served

i j

by the facility in location j X is the array Ixij]

L is the subset of locations of open facilities, to be chosen among all subsets of feasible facility locations

Cij is the cost paid to serve a customer in location i by the facility in location j; it includes transport costs and possible operating costs

a is the fixed cost paid to establish a facility in j

location j

The goal of problem ( 1 ) - ( 3 ) is therefore to find a spatial

arrangement of facilities, L, and an assignment of customers to them, X, which minimize total cost.

A well known property of ( 1 ) - ( 3 ) is that customers are assigned to the least-cost facility. Indeed, for L fixed, problem ( 1 ) - ( 3 ) is separable for each customer location, and the subproblem associated with a customer location i is

min Z xij X jEL i' j

with the trivial solution

The assignment rule ( 4 ) can be given two different interpretations.

If the assignment of customers to facilities is controlled by the decision maker, it is the optimal solution to the total cost

minimizing problem. If customers are free to choose the facility they want, it is a model for rational choice behavior, stating that customers choose according to the minimum-cost criterion.

The last interpretation is the point of departure for the gener- alization which will be developed.

Let it be assumed that facilities belong to the second

class discussed above, so that customers choose facilities, have to travel in order to get served, and pay for transport and

possible operating costs. The behavior implied by (4) is not only rational but also deterministic, since no uncertainty, lack of information, or variation in tastes and preferences is ac- counted for. On the other hand, everyday experience suggests that cost is not the only criterion determining customer choice.

If the assumption of rational choice behavior is kept, it can be said that customers maximize their utility, which is basically a function of costs but also includes many other variables.

Moreover, not all the components of this utility can be measured easily, since they vary greatly among individuals, and possibly among different points in time for the same individual.

A simple model which accounts for the above requirements is as follows. Let it be assumed that the utility of choosing a

facility in location j, for a customer in location i, can be split

in two parts, a measured (or deterministic) one and a random one

u w _{i j} = vij +

Fj

₍₅₎

where

..

u

i j is the total utility of choosing a facility in location j, for a customer in location i

v ij is the measured part of utility

Fj

is the random part of utility

The

9

are random variables with a joint distribution function, j

from which all customers are assumed to draw.

Each customer in i will choose in order to maximize his utility, that is, he will solve the problem

max ii jEL ^ij

But the quantity (6) is a random variable. The expected utility for a customer in i is therefore

Ui (L) = E max G (,EL ij)

where E denotes the expectation operator. The random-utility counterpart of problem ( 1 )

-

⁽³⁾ is therefore

max CP.U.(L)

-

^{C a}

L i ^{1 1} j EL j

Notice that the rnax operator has replaced the m i n operator, since utility maximizing rather than cost minimizing, is used.

Notice also that constraints (2) and (3) have been dropped, and the array X has disappeared, since the assignment subproblem has already been solved by introducing (6) and (7). Problem (8) looks therefore simpler than problem (1)-(3). What price is paid for this "simplicity"? Basically, the linear-integer

programming features of (1 )

-

⁽³⁾ are lost. While problem (1 )

-

⁽³⁾

can be easily solved by dual ascent and related methods (Bilde and Krarup, 1977; Erlenkotter, 1978; Wolsey, 1980), problem (8) cannot. It is a general combinatorial optimization problem for which satisfactory exact algorithms are not known, except for problems of small size. However, the random-utility assumption will be shown to give problem (8) some special properties, which can be exploited to develop f.airly good heuristic approaches.

3. BASIC PROPERTIES OF THE ADDITIVE RANDOM-UTILITY MODEL

In order to operationalize the model loosely introduced in Section 2, some basic notions from random-utility theory are needed. Most results presented in this section can be found in the recent literature (Domencich and McFadden, 1975; Williams,

1977; Daly, 1978; Ben-Akiva and Lerman, 1978).

Some of them are new, and specially tailored to give insight into problem (8)

.

The following notation will be used. An upper case X,Y, etc., is a vector; a lower case x,y, etc., is either a scalar or a vector with all elements equal to x,y, etc. Therefore

Functions of vectors are defined element by element, e.g.,

2 i s a random v a r i a b l e , u s u a l l y t h e j t h e l e m e n t o f a random j

v e c t o r . 2 i s a random v e c t o r , whose d i s t r i b u t i o n f u n c t i o n i s

A l l t h e d i s t r i b u t i o n f u n c t i o n s c o n s i d e r e d i n t h i s p a p e r a r e assumed t o be c o n t i n u o u s and t o have d e r i v a t i v e s o f any o r d e r . The c o n d i t i o n a l d e n s i t y of 2 i s d e n o t e d by F . ( X ) and d e f i n e d

j I

by e q u a t i o n s

Fj _{( X ) =} aF ( X )

ax j

The extreme-value d i s t r i b u t i o n of X

-

i s t h e d i s t r i b u t i o n f u n c t i o n of t h e random v a r i a b l e

max 2 j j

I t i s w e l l known from t h e t h e o r y o f extreme o r d e r s t a t i s t i c s (see Galambos, 1978, f o r i n s t a n c e ) t h a t t h e f o l l o w i n g e q u a t i o n h o l d s t r u e

t h e r e f o r e , from (13) and ( 9 ) i t f o l l o w s

The extreme-value p r o b a b i z i t y f o r element 2 i s t h e p r o b a b i l i t y j

From equation (10) it is easily derived that

Let now the additive random-utility model be introduced.

In order to simplify notation, the subscript labeling customer location will be dropped, since it will be kept constant. From

( 5 ) , the total utility when alternative j is chosen is

Let F(Y) be the distribution function of Y =

{?;I.

Then the

J

distribution function of

6

⁼

^{{ G .}

^is

3

where V = {v.) is the vector of measured utilities.

3

The extreme value distribution for

5

is, according to (14)

and its first moment is

Equation (19) gives the expected utility for a rational customer.

The extreme value probability for alternative j is the

probability that the facility in j is chosen. According to (16) it is given by

where F ₃

.

^(X) is defined by (1 1 )

.

The following propositions state some noteworthy properties of the additive random-utility model.

PROPOSITION 1. (Translation of expected utility)

P r o o f

- OD

-

^{ydF (y}

-

^V) + ^a _{dF (y}

-

^{V) = $ (V)}

+

^a

(The transformation y = x

-

a has been used.)

PROPOSITION 2. (Translational invariance of choice probabilities)

P r o o f

(The transformation y = x

-

a has been used.)

PROPOSITION 3. (Hotelling consistency of expected utility)

P r o o f

[Equations (11) and (20) have been used.]

PROPOSITION 4. (Logit-like representation of choice probabilities)

T h e F u n c t i o n

'4 (W) = exp4 (log W)

is linear homogeneous, a n d

w h e r e

P r o o f :

$ (aw) = exp4 ( l o g W

+

l o g a ) = exp [ $ I ( l o g W)

+

l o g a ] = a Q ( W ) hence l i n e a r homogeneity f o l l o w s

[ P r o p e r t y ( 2 1 ) h a s been used. I

To p r o v e ( 2 5 ) , l e t ( 2 4 ) b e used t o w r i t e t h e e x p e c t e d u t i l i t y i n t h e form

4 ( V ) = l o g $ (eV)

t h e n

[ P r o p e r t y ( 2 3 ) h a s been used. 1

PROPOSITION 5. ( E q u i v a l e n c e w i t h an e n t r o p y maximizing problem)

.4 j

@ ( V ) = max - L qj l o g ^f : L qj = 1

Q j j j

where

Q =

a n d $ . ^{( 0} ) i s d e f i n e d b y ( 2 6 ) I

P r o o f :

Define

H(Q) = - Z q . log 1 j I

j

Then for any Y = {yj }

1 yj Yk Hjk (Q) = -E 2

-

^{< 0}

jk j qj

hence H ( Q ) is concave.

The solution to the mathematical program

satisfies the conditions

qj ⁼ cif j

where v is a Lagrange multiplier and a = e ^-(1

+

^v)

Elimination of a by the constraint yields

and the concavity of H(Q) ensures that this solution globally maximizes H (Q)

.

But from (28)

[The linear homogeneity of $ ( a ) has been used.]

Hence the q obtained by solving the above mathematical program are the same as the ones given by j ( 2 5 ) . Sub- stitution into H(Q) yields

H(Q) = - z q . log _{= log lb(e}

v

_{= Q(V)}

j ^{J $} (eV)

[Equation (24) has been used. ]

PROPOSITION 6. (Nondecreasing submodularity of expected utility)

T h e s e t f u n c t i o n

U(L) = E(max _{i i . )} jEL J

is s u b m o d u Z a r nondecreasing.

P r o o f

The f u n c t i o n

G ( L )

^{= max ii}

j EL j

i s submodular n o n d e c r e a s i n g (Nemhauser, Wolsey, and F i s h e r , 1978)

T h i s p r o p e r t y c a n be s t a t e d a s

G ( s

^U i ^{j ) )}

- ^fi(s) -

^>

^{G ( T}

~ { j ) )

- ^{G ( T )}

^{1 0}

f o r a l l S , T , j

,

^{S C T} -

,

^{j F T}

Applying t h e e x p e c t a t i o n o p e r a t o r t o b o t h s i d e s o f t h e above i n e q u a l i t y o n e g e t s

h e n c e U ( L ) i s submodular n o n d e c r e a s i n g .

B e s i d e s t h e above p r o p o s i t i o n s a n o t h e r p r o p e r t y w i l l b e u s e f u l , s t a t i n g t h e r e l a t i o n s h i p between ( 1 9 ) and ( 2 9 ) . The r e l a t i o n s h i p i s

u ( L ) = l i m $ (V) v~

+

^{- w}

1FL

Most p r o p o s i t i o n s above a r e s e l f - e x p l a n a t o r y and need o n l y a few comments.

Propositions 1 and 2 state that choice behavior is unaffected by shifts in utilities. In other words, the "zero" of the utility

scale may be set arbitrarily.

Proposition 3 is perhaps the most important one, since it states the integrability condition for choice probabilities. The importance of this property is well known, and it has been dis- cussed by many authors (for instance, Williams, 1977; and Daly,

1978).

Proposition 4 gives a useful representation of choice

probabilities, and shows how the general additive random-utility model is related to the well-known logit model (included as a special case). Equation (25) qeneralizes a result first obtained by McFadden (1978) under much more restrictive assumptions.

Propositions 5 and 6 state some useful properties to analyze the optimal location problem. Equation (27) relates the additive random-utility model to the wide and well-known class of entropy maximizing models (Wilson, 1970; Wilson and Senior, 1974;

Willekens, Por, and Raquillet, 1979). The submodularity property relates the location problem discussed in this paper with the problems analyzed by Nemhauser, Wolsey, and Fisher (1978), for which some approximate optimization results have been produced.

4. BASIC PROPERTIES OF SUBMODULAR SET FUNCTIONS

Most results given in the preceding section were aimed at providing theoretical insight and economic justifications for the behavioral models based on the additive random-utility

assumption. At least one of them (namely proposition 6), however, is specially tailored on the more operational goal of providing techniques to solve problems of type (8). Since the submodularity property is the key premise to all subsequent results (the non- decreasing property is not essential), it will be shown to hold true for functions built like (8). First, since the functions Ui(L) are submodular because of proposition 6 and since any

positive linear combination of submodular functions is sub- modular (Nemhauser, Wolsey, and Fisher, 1978), the function

is submodular. Let the set function

be defined. Then

Because of submodularity of g(L), the inequality

holds true for all L,T,j, L

-

^{C T, j} T. Subtracting a on j both sides of (34) and using (33), the inequality

for all L,T,j, L

-

^{C T, j}

^P

^T follows, thus stating submodularity of G(L). The general optimization problem is therefore to find the unconstrained maximum of a submodular function

max = G(L) L

If G(L) were also nondecreasing the solution to problem (36) would be trivial, i.e., a facility should be established in a22

possible locations. However, the presence of (usually nonzero) fixed charges a prevents the nondecreasing property from being true. By means of a j economies of scale are introduced, and

j

'

the higher their value the smaller will be the number of facili- ties to be established.

The properties of submodular set functions, in the light of cornbinatorial optimization, have been thoroughly analyzed recently by Cornuejols, Fisher, and Nemhauser (1977); Nernhauser, Wolsey, and Fisher (1978); Fisher, Nemhauser, and Wolsey (1978);

and Wolsey (1980). The reader is referred to their work for further results. Here only a few statements will be needed, which do not go far beyond alternative definitions of the sub- modularity property itself.

Let the following quantities be defined as

pj (S) = G(SUfj1)

-

^G(S) the incremental value of

adding element j to the (37) set S, for j

9

^S

The p.(S) give information on how the value of the objective function changes when the solution set S is augmented by one 3

location j.

Often use will be made of the change in the value of the objective function due to the deletion of one location j E S .

From equation (37) it follows that this change is given by

The statements on submodular set functions which will be used in the rest of the paper are summarized in the following proposition.

PROPOSITION 7.

E a c h o f t h e f o l l o w i n g i n e q u a l i t i e s d e f i n e s a s u b m o d u l a r s e t f u n c t i o n

pj(S)

2

^p.(T) ₃

,

f o r a l l S

-

^{C T a n d j} 9 ^{T (38}

G(T) < G ( S ) + ^{1 pj(S)}

-

^Z pj(SUT-{j])

j ET-S jES-T

/

⁽³⁹⁾

f o r aZZ S and T

f o r aZZ S

-

C T ⁽⁴⁰⁾

G(T) <G(S)

-

^{1 p.(S-{j))} f o r aZZ T

-

^C

^s

jES-T 3

The proof of the above propodition is found in the four refer- ences listed above.

For an intuitive feeling of how a submodular set function looks, compared to the usual functions of sets of real numbers, one could note that the submodularity property is a kind of generalization of concavity. Indeed, if the quantities p.(S)

3 are considered as analogous to derivatives, property ( 3 8 )

suggests the notion of decreasing return to scale; i.e., adding one more element to the solution set is always more beneficial for smaller sets than for bigger ones.

In order to devise techniques to solve problem (36), a first step is to impose some restrictive requirements on possi- ble solution sets. One is thus led to look for some analog of

ZocaZ m a x i m a , i.e., sets which are optimal within some suitably

defined neighborhood.

The simplest local properties for maxima are stated in the following definitions

DEFINITION 1.

A weak ( o r f i r s t o r d e r ) ZocaZ maximum o f a s e t f u n c t i o n G(S) i s a s e t S s u c h t h a t

I

pj(S)

5

⁰

,

f o r aZZ j E

s

(

~ ~ ( s - { j ] ) l O

,

f o r a t 2 j E S

DEFINITION 2.

A s t r o n g f o r s e c o n d o r d e r ) Zocat maximum o f a s e t f u n c t i o n G ( S ) i s a s e t S s u c h t h a t

pj(s) = min pk(sU{jl-{kl) 2 0

,

f o r a t t j P S kESUC j

l

2

o r , e q u i v a t e n t t y

pj(S- {jl) = max %(S-{jl) 1 0

,

f o r a t t j E S kqs-C j

1

Stated in words, a weak local maximum is a set which cannot be improved by the addition or deletion of a single element. A

strong local rnaxir-um is a weak local maximum too; moreover, it cannot be improved by any paired interchange between two elements.

Of course generalizations to higher order local maxima can be devised, requiring stability under interchange of n elements, for n > 2. But the numerical results shown later seem to suggest that going beyond the second order is seldom needed.

If G(S) is submodular, the following propositions can be stated for weak local maxima.

PROPOSITION 8. (Dominance over all supersets) I f S i s a weak ZocaZ maximum and S

- ^c

^{T, t h e n}

Proof

from (42) p.(S) < 0, for j E T - S ; substitution into

3

-

(40) yields

PROPOSITION 9. (Dominance over all subsets) I f S i s a weak l o c a l maximum and T - ^C S, t h e n

Proof

from (42) pj(S-{jI) 2 ⁰

,

^{for j E S - T ;}

substitution into (41) yields

The above two simple propositions surprisingly widen the neighborhood dominated by a weak local maximum. What they state

is that if a weak local maximum is detected, then any of its subsets or supersets can be dropped from further search.

Two further propositions will be useful to build rules of improvement for a possible trial solution.

PROPOSITION 10. (Nondecreasing point detection)

I f pj(~) 2 ⁰ f o r some j

9

^T, t h e n no S c T is a weak l o c a Z maximum

Proof

from (38)

,

^Pj(s)

2

^Pj(T)

-

⁰

,

^{for some j}

9

^S

which contradicts (42).

PROPOSITION 1 1 . (Nonincreasing point detection)

~f p j ( ~ - {j}) 5 ⁰ f o r some j E T, t h e n no S

-

3 T _{i s} a weak l o c a l maximum

Proof

since S - {j)

-

3 T - {j), from (38)

pj (S

-

^{{j)) pj (T- {jl)}

⁵

^{0, for some j E S,}

which contradicts (42)

.

Propositions 10 and 1 1 are useful for building possible tree- search algorithms. If the search goes downward, that is

building small and smaller subsets, then Proposition 10 may be used as a rule to stop searching in a subset. If the search goes upward, that is building bigger and bigger supersets, then Proposition 1 1 may be used as a rule to stop searching in a superset.

Proposition 8-1 1 refer to weak local maxima. One would hope to find some more restrictive conditions when strong local maxima are used. Unfortunately, this does not seem to be the case. No further property of strong local maxima has been found as yet, except those stated in the definition.

5. SOME ALGORITHMS

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 . a}

i j I j

be defined over all real nonnegative vectors X = x , where the

4 . (V) ₁ are the functions defined by equation (1 9)

.

^{Then, if}

it 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 ) and

jgs

^I

p 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 p}

i f p i ( S ) >

-

p ( s - { k } ) r e p l a c e S by s U { i } k

i 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

+

y

a 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) where

Comparison 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

- _a

Since, 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 a}

i 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 a}

L 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 T a b l e 1 . The r e s u l t s f o r H e u r i s t i c s 1 and 2 have been produced w i t h t h e s t a n d a r d s t a r t ( 5 0 ) .

Table 1 i s s u r p r i s i n g i n many ways. F i r s t o f a l l , i t shows t h e u n e x p e c t e d power o f H e u r i s t i c 1 . I f t h e 2nd and 4 t h column a r e compared, it i s s e e n t h a t H e u r i s t i c 1 f a i l e d t o f i n d t h e e x a c t s o l u ' t i o n o n l y f o r a f i x e d c h a r g e of 2,500. Even i n t h a t c a s e , 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 (boxed i n T a b l e 1 ) i s v e r y c l o s e t o t h e o p t i m a l one. Comparison o f t h e 3rd and 4 t h column i s a l s o r e v e a l i n g . The 3rd column shows t h e r e s u l t s

o b t a i n e d i n t h e f i r s t s t a g e o f t h e INTLOC a l g o r i t h m , t h a t i s t h e b e s t lower bounds produced b e f o r e e n t e r i n g t h e branch-and-bound r o u t i n e . Except f o r a f i x e d c h a r g e v a l u e o f 500 ( f o r which a l l methods g i v e t h e o p t i m a l s o l u t i o n s ) , t h e s e v a l u e s a r e a l w a y s

nonoptimal. Even f o r a f i x e d c h a r g e o f 2,500 t h e v a l u e found w i t h H e u r i s t i c 1 , a l t h o u g h n o n o p t i m a l , i s c l o s e r t o t h e o p t i m a l o n e t h a n t h e one produced w i t h t h e f i r s t s t a g e o f INTLOC.

H e u r i s t i c 1 seems t h e r e f o r e t o d e f i n i t e l y o u t p e r f o r m t h e f i r s t s t a g e of INTLOC.

H e u r i s t i c 1 h a s a l s o been t r i e d w i t h s t a r t s d i f f e r e n t from ( 5 0 ) , and t h e r e s u l t s ( n o t r e p o r t e d h e r e ) have n o t a l w a y s been s o good. The p r o c e d u r e o f t e n t e r m i n a t e d on n o n o p t i m a l l o c a l maxima. However, even i n t h e w o r s t c a s e s , t h e comparison w i t h t h e f i r s t s t a g e o f INTLOC h a s always been i n f a v o r o f H e u r i s t i c 1 .

The second i m p o r t a n t f a c t shown by T a b l e 1 i s t h e e f f e c t i v e - n e s s of H e u r i s t i c 2 , which seems t o work much b e t t e r t h a n e x a c t a l g o r i t h m s . T h i s e f f e c t i v e n e s s i s a c t u a l l y much h i g h e r t h a n shown i n t h e t a b l e . H e u r i s t i c 2 h a s been s y s t e m a t i c a l l y t e s t e d w i t h many d i f f e r e n t s t a r t s , and it n e v e r f a i l e d t o r e a c h t h e optimum f o r e v e r y f i x e d c h a r g e v a l u e .

The comparison of computing times i s a l s o r e v e a l i n g . F o r t h e e x a c t s e a r c h methods, o n l y t h e CPU t i m e s f o r t h e method b a s e d on s u b m o d u l a r i t y bounds h a v e been r e c o r d e d , b u t t h e y a r e o f t h e same o r d e r of magnitude a s f o r INTLOC. The a s t e r i s k s i n t h e 6 t h column i n d i c a t e t h a t t h e t r e e s e a r c h w a s t o o s l o w t o r e a s o n a b l y w a i t f o r a s o l u t i o n , and t h e a l r e a d y a v a i l a b l e s o l u t i o 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?

8. CONLCUDING COMMENTS AND ISSUES FOR FUTURE RESEARCH

The main aim of the last sections of this paper has been to show the effectiveness of simple ascent heuristics when applied to nonconventional facility location problems (i.e., whose

objective function is based on random-utility theory!. The key property leading to these results seems to be the submodularity of the objective function. The problem formulation itself may lack some realistic features, since actual location problems

often have constraints not considered here, like budget and size constraints. However, the successful solution of the simple problem is the main step towards solving more complex ones.

Related work has shown (Leonardi and Bertuglia, 1981) that effective ascent heuristics can be built for problems with capacity constraints as well. Future progress can be expected

in solving problems of the type recently explored by Coelho (1980b), including fixed charges, capacity constraints and p-median type

constraints (that is, constraints on the number of facilities to be established).

Another strand of future research is the exploration of different objective functions. The numerical results reported in this paper are based on the special assumptions introduced in Section 6. Other forms of random-utility distributions could be tried, leading to other objective functions. However, the submodularity property would still hold because of Proposition 6.

Further mathematical investigation is also required, as already stated in Section 7. The reason why the proposed heuristics outperform other methods is far from being fully understood, although some theoretical results for related

simpler problems (Wolsey, 1980) seem to suggest that such a good performance is not too surprising. This paper is only an

exploratory one, showing a way which might be worth following.

APPENDIX: INPUT DATA USED FOR NUMERICAL TESTING ( T u r i n High S c h o o l s )

Travel times (minutes) on public transport

High s c h o o l t o d i s t r i c t demand

D i s t r i c t 1 2 3 4 5 6 7 8 9 1 0 1 1 1 2 13 14 15 1 6 17 1 8 1 9 2 0 2 1 2 2 23

(No. s t u d e n t s ) from d i s t r i c t

Travel time discount rate: a = 0.194.

REFERENCES

Beaumont, J.R. (1979) Some I s s u e s i n t h e A p p l i c a t i o n o f Mathemat-

i c a l Programmi~zg i n Human Geography. WP-256. Leeds, UK:

School of Geography, University of Leeds.

Beaumont, J.R. (1980) Spatial interaction models and the location- allocation problem. J o u r n a l o f R e g i o n a l S c i e n c e 20:37-50.

Ben-Akiva, M., and S. R. Lerman (1978) Disaggregated travel and mobility-choice models and measures of accessibility. Pages

654-679 in S p a t i a l I n t e r a c t i o n T h e o r y and P l a n n i n g M o d e l s , edited by A. Karlqvist, L. Lundqvist, F. Snickars, and

J.W. Weibull. Amsterdam: North Holland Publishing Company.

Bilde, O., and J. Krarup (1977) Sharp lower bounds and efficient algorithms for the simple plant location problem. A n n a l s o f D i s c r e t e M a t h e m a t i c s 1:79-97.

Coelho, J.D. (1977) The Use o f M a t h e m a t i c a l O p t i m i z a t i o n Methods i n Model Based Land Use P l a n n i n g : An A p p l i c a t i o n t o t h e New Town t o S a n t o Andre. Ph.D. Thesis. Leeds, UK: University of Leeds.

Coelho, J.D. (1979) A locational surplus maximisation model of land use plan design. Pages 48-60 in London Papers i n R e g i o n a l S c i e n c e : D e v e l o p m e n t s i n Urban and R e g i o n a l A n a l y s i s , edited by M . J. Breheny. London: Pion.

Coelho, J.D. (1 980a) ~ p t i r n i z a c ; i o , ~ n t e r a c c z o E s p a c i a l e ~ " e o r i a d o Comportamento ( O p t i m i z a t i o n , S p a t i a l I n t e r a c t . i o n , and B e h a v i o r T h e o r y ) . Nota No. 5 .

de Estatistica e Aplicac~es, Dept. Matematica Aplicada, Faculdade de Ciencias de Lisboa.

Coelho, J.D. (1980b) A L o c a t i o n a l S u r p l u s M a x i m i z a t i o n A p p r o a c h t o P u b l i c F a c i l i t y . L o c a t i o n . Nota No. 15. Lisbon,

Portugal: Centro de Estatistica e ~plicacges, Dept.

Matematica Aplicada, Faculdade de Ciencias de Lisboa.

Coelho, J.D., and H.C.W.L. Williams (1978) On the design of land use plans through locational surplus maximization. P a p e r s o f t h e R e g i o n a l S c i e n c e A s s o c i a t i o n 40:71-85.

Coelho, J.D., and A.G. Wilson (1976) The optimum location and size of shopping centres. R e g i o n a l S t u d i e s 4:413-421.

Cornuejols, G., M.L. Fisher, and G.L, Nemhauser (1977) Location of bank accounts to optimize float: An analytic study of exact and approximate algorithms. Management S c i e n c e 23:789-810, Daly, A.J. (1978) Some developments in transport demand modeling.

Pages 319-333 in S p a t i a l I n t e r a c t i o n T h e o r y a n d P l a n n i n g M o d e l s , edited by A. Karlqvist, L. Lundqvist, F. Snickars, and J.W. Weibull. Amsterdam: North Holland Publishing Company.

Domencich, T., and D. McFadden (1975) Urban T r a v e l Demand: A B e h a v i o r a l A n a l y s i s . Amsterdam: North Holland ~ublshing Company.

Efroymson, M.A., and T.L. Ray (1966) A branch-bound algorithm for plant location. O p e r a t i o n s R e s e a r c h 14:361-368.

Erlenkotter, D. (1978) A dual-based procedure for uncapacitated facility location. O p e r a t i o n s R e s e a r c h 26:992-1009.

Erlenkotter, D., and G. Leonardi, A l g o r i t h m s f o r S p a t i a l

I n t e r a c t i o n - B a s e d F a c i l i t y L o c a t i o n M o d e l s . Forthcoming Working Paper. Laxenburg, ~ustria: International Institute

for Applied Systems Analysis. (Presented at the IIASA Norm- ative Location Modeling Task Force Meeting, June 10-13, 1980.) Ermoliev, Y., G. Leonardi, and J. Vira (1981) T h e S t o c h a s t i c

Q u a s i - G r a d i e n t Method A p p l i e d t o a ~ a c i l i t y L o c a t i o n P r o b l e m . WP-81-14. Laxenburg, Austria: International Institute for Applied Systems Analysis.

Fisher, M.L., G.L. Nemhauser, and L.A. Wolsey (1978) An analysis of approximations for maximizing submodular set functions-11.

M a t h e m a t i c a l Programming S t u d i e s 5 : 73-87.

Galambos, J. (1978) T h e A s y m p t o t i c T h e o r y o f E x t r e m e O r d e r S t a t i s - t i c s . New York: John Wiley.

Gumbel, E.J. (1958) S t a t i s t i c s o f E x t r e m e s . New York: Columbia University Press.

Leonardi, G. (1975) Un nuovo algoritmo per il problema della localizzazione ottimale dei servizi urbani (A new algorithm for the optimal location or urban services). P r o c e e d i n g s o f t h e C o n g r e s s o f t h e I t a l i a n O p e r a t i o n s R e s e a r c h A s s o c i - a t i o n , Milan.

Leonardi, G. (1978) Optimum facility location by accessibility maximizing. E n v i r o n m e n t and P l a n n i n g A 11:1287-1305.

Leonardi, G. (1980a) A U n i f y i n g Frameuork f o r P u b l i c F a c i l i t y L o c a t i o n ProbZems. WP-80-79. Laxenburg, Austria:

International Institute for Applied Systems Analysis.

(Forthcoming in E n v i r o n m e n t and P l a n n i n g A . )

Leonardi, G. (1980b) A M u l t i a c t i v i t y L o c a t i o n Model w i t h

A c c e s s i b i Z i t y - and C o n g e s t i o n - S e n s i t i v e Demand. WP-80-124.

Laxenburg, Austria: International Institute for Applied Systems Analysis.

Leonardi, G., and C.S. Bertuglia (1981) O p t i m a l High S c h o o l L o c a t i o n : F i r s t R e s u l t s f o r T u r i n , I t a l y . WP-81-05.

Laxenburg, Austria: International Institute for Applied Systems Analysis.

McFadden, D. (1978) Modelling the choice of residential location.

Pages 75-96 in S p a t i a l I n t e r a c t i o n T h e o r y and P l a n n i n g

M o d e l s , edited by A. Karlqvist, L. Lundqvist, F. Snickars, and J.W. Weibull. Amsterdam: North Holland Publishing Company.

MacGill, S.M., and A.G. Wilson (1979) Equivalences and similar- ities between some alternative urban and regional models.

S i s t e m i U r b a n i 1:l-40.

Nemhauser, G.L., L.A. Wolsey, and M.L. Fisher (1978) A n analysis of approximations for maximizing submodular set functions-I.

MathematicaZ Programming 14:265-294.

Provincia di Torino ( 1 978) La PopoZazione S c o l a s t i c a de ZZe S c u o l e Medie S u p e r i o r i S t a t a Z i i n P r o v i n c i a d i T o r i n o . (The

Student Population in the Public High School in the Province of Turin. ) Turin, Italy: Provincia di Torino.

Spielberg, K. (1969) Algorithms for the simple plant-location problem with some side conditions. O p e r a t i o n s R e s e a r c h

17:85-111.

Van Lierop, W., and P. Nijkamp (1979) A utility framework for interaction models. S i s t e n i U r b a n i 1:41-64.

Van Roy, T.J., and D. Erlenkotter (1980) A Dual-Based ProceEure f o r Dynamic F a c i Z i t y L o c a t i o n . WP-80-31. Laxenburg, Austria: International Institute for Applied Systems Analysis.