Dual search procedures for the exact formulation of the simple plant location problem with spatial interaction

Pergamon

LocnfionScience,Vol.4, No. l/Z, PP.83-100,1996 Copyright0 1996Else&r ScienceLtd Printedin Great Britain.All rights reserved 0966-8349196 $15.00+ 0.00

PII: so%6-8349(%)ooo13-7

DUAL SEARCH PROCEDURES FOR THE EXACT FORMULATION OF THE SIMPLE PLANT LOCATION PROBLEM WITH SPATIAL INTERACTION KAJ HOLMBERG LinkCiping Institute of Technology, Department of Mathematics, S-581 83 LinkGping, Sweden

KURT JiiRNSTEN Norwegian School of Economics and Business Administration, Institute of Finance and Management Science, Helleveien 30, N-5035 Bergen-Sandviken, Norway (Received for publication April 1996)

Abstract-We present an exact formulation of the simple plant location problem with spatial interaction. The formulation is based on the derivation of the spatially interactive travel behavior, often described. through a spatial “gravity” model. By using techniques from separable programming, we can derive an exact formulation of the problem, i.e. without using any approximation of the function that describes the spatially interactive travel behavior. The resulting model is a pure zero-one location model. We also present a dual ascent and adjustment procedure for the linear programming relaxation of the exact formulation of the simple plant location problem with spatial interaction. Furthermore, we present a solution method based on Lagrangean relaxation and subgradient optimization. The two dual search procedures are computationally tested on problems with a maximal size of 50 x 100. Copyright 8 1996 Elsevier Science Ltd Keywords: Location, spatial interaction, mixed integer programming, dual ascent.

1. INTRODUCTION Research on facility location models has mainly dealt with location models in which the allocation rule is to assign demand to its nearest (or more precisely, its minimum cost) facility. This assignment rule is certainly appropriate when the cost of both establishing the facilities and the transportation of goods are paid by the producer. However, for many public facility location problems, the clients are free to make their own choice of facility and it can be observed that the cost minimization based allocation rules do not give a good description of real life behavior in these situations. In order to get a model that is better suited for these situations, many researchers have formulated facility location models in which the spatially interactive travel behavior is described through a spatial “gravity” model. Location models including spatially interactive travel behavior have been presented by Coelho and Wilson (1976), Leonardi (1978), Beaumont (1980), Erlenkotter and Leonardi (1985) and Jacobsen (1986). The foundation for this behavior in location models has been shown by the use of random utility theory by Leonardi (1983). The resulting models are described as pure location models including an entropy term. The models are thus nonlinear integer 83

84

K. HOLMBERGand K JGRNSTEN

programming models. Solution methods for this class of problems have been presented in some of the papers above, and in Coelho (1980) and Jacobsen (1988). In this paper, we will show that by using the classical way of deriving the gravity model without doing any approximation, we obtain a nonlinear integer programming model, which can be reformulated into a linear binav programming problem, closely related to the simple plant location problem. We also describe a dual ascent procedure for solving the problem, which increases the dual variables in small steps, while keeping under control, via the complementary slackness conditions and dual feasibility, the (monotone) advance towards the optimum, i.e. towards primal feasibility. We then use a dual adjustment procedure, which decreases dual variables that block further improvement of the dual ascent procedure, and then attempt to increase other dual variables. The whole method is a heuristic, but allowing for the simultaneous change of several dual variables, one could guarantee covergence towards the optimal objective function value of the LP-relaxation. Another possibility is to use a dual subgrudient optimization method, based on Lagrangean relaxation, which has asymptotic (not monotone) convergence towards the optimal objective function value of the LP-relaxation. By computational tests, we compare the practical performance of the dual ascent and adjustment method with dual subgradient optimization, and we find that they complement each other in an interesting way. In case there is a remaining duality gap between the objective function value of the LPrelaxation and the optimal objective function value of the integer problem (occurring when the LP-relaxation returns to a fractional solution) after the dual search procedure has converged, one can use branch-and-bound to close it. In section 2, we derive the simple plant location problem with spatial interaction. Section 3 describes the linearization procedure and gives the final exact formulation. In section 4, Lagrangean relaxation and subgradient optimization are described. In section 5, the LPrelaxation of the plant location problem with spatial interaction is given, together with its LP-dual. In section 6, we show how to solve (optimally) the remaining parts of the dual, when certain dual variables are fixed. We also use the complementary slackness conditions and the structure of the problem to find a practical characterization of the optimal solution and to get a useful optimality criterion. We describe the dual ascent procedure, based on the same principles as that of Erlenkotter (1978), discuss the dual adjustment part and finally give the complete algorithm. Section 7 describes our computational tests on random test problems of maximum size 50 x 100 and section 8 contains some conclusions. The contribution of this paper is a new location model that can be efficiently solved. The model has the advantages of older models, both when it comes to generality and applicability of the model and possible solution methods. We have developed a new dual ascent method, similar in spirit to that of Erlenkotter (1978), but technically quite different, as the model is different. We have also applied Lagrangean relaxation and subgradient optimization. Computational tests show that with these two methods available, it is possible to solve problems larger than previously have been solved in the literature in a reasonable time. 2. THE SIMPLE PLANT LOCATION PROBLEM WITH SPATIAL INTERACTION Suppose that we have m possible locations and n demand points. If plant i is opened, a fixed cost, ui, occurs. At demand point j the demand is We Trips will emanate from the demand points and arrive at the opened plants so that the demand is satisfied. The

Dual search procedures

85

transportation costs for one trip from demand point j to plant i are c+ The decision variables are +=number

of trips from demand point j to plant i, i=l,. . . ,m, j=l,. . . ,n. 2i=

1 if plant i is opened

i=l,...,m.

0 if not

In the simple plant location problem, SPLP, the objective is to minimize the total cost, which is the sum of all fixed costs and all transportation costs. The SPLP has a linear programming relaxation that often yields integer solutions. This is used in some existing efficient solution methods; see, for example, Erlenkotter, 1978, Bilde and Krarup, 1977 and Korkel, 1989. If we identify every single trip in a solution, we obtain a set of microstates, while the aggregated solutions x are the macrostates. If the trips are made by individuals, one cannot expect that the solution that occurs is the one with the minimal total cost. In situations in which the clients are free to make their own choice of facility, there exists empirical evidence that not all clients will select their nearest facility. In order to model this spatially interactive behavior, we follow Wilson’s (1967) approach, to find the macrostate that is associated with the maximum number of microstates. It can be argued that the solution that maximizes the number of microstates is the most probable solution. Including this maximization in the objective function is equivalent to m

n

min JJ JJ

m

Xij

or min 2

n

C In(x, !).

Our “exact” formulation of the simple plant location problem with spatial interaction, abbreviated SPLPS, is the following. Identify ‘u*,where

S.t.

f

Xij=Wj

j=l,...,

n

(1)

i=l

X,i-WjZi
i=l,...,

Xij2 0, integer zi=O/l

Vi

mJ=l,..., Vi,j

n

(3 (3) (4)

(Pl)

86

K. HOLMBERG and K. J6RNSTEN

where y reflects the effect of the spatial interaction relative to the transportation and location costs. Constraint set 1 ensures that the demand will be met, while constraint set 2 prohibits trips to unopened plants. Note that instead of using a model with budget constraints as is done by Erlenkotter and Leonardi, we have used a model in which we weigh together the two objectives, cost minimization and the maximization of the number of microstates. The reason for this choice of formulation is that it will in most practical applications be of interest to determine how the travel pattern behaves as a function of the total costs. In previous work, Stirling’s approximation has been used to approximate the logarithmic part of the objective function. This yields the simple plant location models with spatial interaction studied by Coelho (1980), Erlenkotter and Leonardi (1985) and Jacobsen (1986), with the same constraints as SPLPS and the following entropy-type objective function.

UC-in

2i

i=l

j-1

(xijln(x,i)-x,i)+y

(2 i i=l

j=l

+++i=lf

aizi >

This yields a mixed integer nonlinear programming problem. In our exact formulation, SPLPS, we do not use the approximation above, instead we keep the exact objective function. Thus, no error is introduced at this point, contrary to what is usual. It should also be observed that, in the derivation, all arguments used are purely combinatorial and the objective function is defined only in integer points. SPLPS is a pure integer problem with a non-linear objective function. In general, such a problem is very difficult to solve. This can account for the use of the approximate formulation presented above. However, due to the structure and separability of the exact nonlinear objective function, we can use another technique to obtain an exact formulation that is a pure binary linear integer problem. We also claim that it might in fact be easier to solve the integer problem the way we propose than to solve the approximate problem as done by Erlenkotter and Leonardi (1985) and Jacobsen (1988).

3. EXACT LINEARIZATION OF THE NONLINEAR OBJECTIVE FUNCTION AND FORMULATION OF THE SIMPLE PLANT LOCATION WITH SPATIAL INTERACTION AS A PURE BINARY PROGRAMMING PROBLEM The reformulation that gives us a pure binary model that is equivalent to the exact simple plant location model with spatial interaction is derived by using techniques from separable programming. The same idea has been used to derive an exact gravity model and is computationally tested in Holmberg and Jornsten (1989). The reformulation is obtained by linearizing the objective function for each variable Xijin the interval O+,
VkgK+

(Letx$-‘)=OVi,j.O!=l and ln(l)=O, so ciil=ln($‘!).)

Dual search procedures

87

Assuming that the break points are equidistant, i.e. xij(k)-~~-l)=djjVk, we can show the following, which indicates convexity of the resulting objective function and is very important for the development of the dual ascent method. TheO?W?l

1: cijk>c@k_l i=l,...,

m,j=l,...,

n, kc2 )...) Wj.

For Vi, vj, VkEKii c~jkk=(ln(~~‘!)-ln(x~-“!))/(~~‘-~~-’))=(ln((~~-‘) +~~)~)-ln(x~~“!))/d~j~((ln(x~~1+d~j)+ln(x~~”+d,,-~)+.~..+ln(2)) -(ln(x~~‘~)+ln(x~~1’-1)+...+ln(2)))/dij=(ln(x~-1~+d,i)+ln(x~-1~+~ij-~) + . ..+

l))/d,, since all the other terms cancel.

ln(xg-‘)+

so for Vi, Vj, Vk&j, k>l,

we have Cijk-Cijk_-l=(lu(X~-l)+d,)+ln(~~-l)+dij-l)

+...+lu(~~-“+l))/d~~-(lu(~,lik~*’+dij)+lu(~~-*’+dii-l)+...+~u(~~-*)

+I))/d~j~(ln(x~~‘~+d~)+ln(x~~“+dij-l)+...+ln(x~~‘~+l)-(ln(x~~2~+~ij) +~n(~~~-*‘+d~~-l)+...+lu(~~-*‘+l))/dij=((lu(~~-‘)+d,,)-lu(~~-*~+dij)) +(~n(X~~‘~+d~-1)-1n(~~~*~+d~~-1))+...+(1~(~~~”+1)-~u(~~~*~+1)))/dij >O, since ln(X~~-“+lij)-ln(x~-*‘+l,)>Otll,=l,...,d,,

due to the simple fact that In(x) is a strictly increasing function. ? ? Now we can peform the substitution

&j=

c

&jkvi,j,

k6,

where x# is the amount of xr’ -xr- ‘) that is used. Let &jk”x$+’-xg- ‘) Vi,j,k.(In the equidistant case we have Uijk=dijVi,j,k.) The linearized problem then becomes: Identity ZL *Lwhere

i=l

S.t.

j=l

2

f i=l

c ‘~1

keK,,

j=l

j=l,...,

xijk’wj

>

i=l

n

(1)

k=l

xiik-wjzi
k=l,...)

Wj,i=l,...,

xijk< uijk Vi,j, k Vi,j, k

xijk2 0 and integer zi=O/l

Vi

m, j=l,...,

n,

(2)

(3) (4) (3

(P2)

K. HOLMBERG and K JORNSTEN

88

The error introduced (i.e. the difference between V* and VT) can only occur if the linearized objective function does not coincide with the exact objective function of SPLPS at some integer point (since non-integral points are not feasible). Clearly, this difference decreases if the distance between the breakpoints, dti, decreases. This difference disappears completely if all integer points are used as break points. Technically, we get xp)=k, Kc= { 1,2,. . . , Wj) Vij, u+=d,=l and c,i,=ln(k!)-ln((k-l)!)=ln(k), k=l,..., wi. The linearization is now exact at every point in which the nonlinear objective function is defined. That is, we have managed to reformulate the simple plant location problem into a pure binary problem of the following form, denoted by SPLPSE. Identify v*, where m

i=l

s.t. f i=l

&jk-zi
2

C Eiikxiik+C yaizi

i-

j=l,...,n

xijk’wj

k=l,...,wj,

Zi=O/l Eijk’cijk

i=l

k-1

(1)

k=l

x,k=o/l

where

j=l

m

w,

”

V*=min C C

i=l,...,

m,

j=l,...,

n,

W)

(2)

Vi,j,k

(3)

Vi

(4)

ycij=ln(k) + YCik

This problem is very similar to the standard simple plant location model. It is also a strong formulation, since the coefficients in the capacity constraints are all one. This is favourable, since it probably means that the gap between the linear programming relaxation and the integer program, which is heavily dependent on the coefficients in the logical constraints, is not too big.

4. LAGRANGEAN RELAXATION WITH SUBGRADIENT OPTIMIZATION Lagrangean relaxation is a well known technique employed for solution of integer and mixed-integer problems; see, for example, Bazaraa and Goode (1979) Geoffrion (1974), Fisher (1981) and Shapiro (1979). Subgradient optimization (see, for example, Shor, 1985) is very often used in conjunction with this approach. The main procedure is to fix the dual variables (Lagrange multipliers) of certain constraints, solve the Lagrangean subproblem (which will be relatively easy), update the values of the dual variables by taking a step in the direction of a subgradient of the dual function and solve the subproblem again. This is repeated iteratively. The subgradient is obtained by inserting the solution obtained from the subproblem into the relaxed constraints and taking the difference between the left-hand sides and the right-hand sides. The step length is determined by a standard step length formula (Held et al., 1974; Poljak, 1969). In each iteration, the subproblem yields a lower bound on v*, but no upper bound is obtained. Convergence can only be expected to be asymptotic. Geoffrion (1974) states that maximizing

Dual search procedures

a9

the Lagrangean relaxation is equivalent to solving the problem with its feasible set replaced by the feasible set of the relaxed constraints intersected with the convex hull of the remaining constraints. Geoffrion (1974) also demonstrated that if the Lagrangean relaxation has the “integrality property” (meaning that the problem yields integer solutions even if the LPrelaxation of the problem is solved), it gives no better lower bound that the LP-relaxation of the whole problem. Relaxing constraint set 1 (and fixing the dual variables Ej to 5.j) yields a subproblem that is separable into one part for each i, containing only one possible location. This means that we replace the demand constraints with linear penalties. The optimal objective function value of the subproblem, h(E), can be obtained as described below.

h(6)=

f

gi(E)+

i=l

i

CrjWj

j=l

where, Vi, n

gi(i)=min

j=l

xijk-Zi
w,

C C

Eiik~iik+y~i~i

k=l

k=l,...,wj,

j=I,...,n

(P4)

xijk=O/I vjik Zi=O/l where Eijk=cijk-&ja (P4), a single facility SPLP, is trivially solvable, by the following considerations. For each i we have two possibilities, Zi=O or ti=l. If ZiZi=Owe get .xijk=OVj,Vk, which yields the objective function value 0. We then assume that Zi=l, and choose the resulting solution only if it is better than 0. Thus, for Vi, we set &(6)=mit&xii,

S.t.Xijk=O/lVj,k

which is trivially solvable, by setting all n’s with a negative objective function coefficient equal to 1 and all the others equal to 0. We get &(c)=min(O,eijk) and

(min(O,&)) +yai >>

with x and z set accordingly. We also note that this subproblem possesses the integrality property, so it cannot produce a bound better than that of the LP-relaxation, vLp. Then we solve the Lagrange dual where uLp < 9 *. This is done by subgradient optimization as described above. The subgradient, 4, which is used as search direction for u, is the difference between the demand and the actual number of trips from the demand point and is obtained as

K. HOLMBERG and K. J6RNSTEN

90

tj=wj_

z i=l

2

fukVj

k=l

where f is the solution of the subproblem. We then update Crby letting c(i=&+ r&jVj, where the steplength r is obtained by the formula

where v is the best upper bound on 2,* known and 0 < E < 1 I 2 -E (E is a constant). 5. THE LP-RELAXATION The linear programming relaxation of the exact formulation, P3, has an optimal objective function value vLp and contains the objective function and constraint set 1 and 2 of P3, with constraint sets 3 and 4 replaced by ziO

Vi

(3)

vi,j,k

(4)

Vi

(5) Let us introduce the dual variables aj for constraint set 1, &k for constraint set 2 and dj for constraint set 3. The LP-dual will then be to identify uLp, where

uLp=max i

set. aj-/$jk
5

2

j=l

k=l

WjCij-

f

j=l

i=l

k=l,...,wj,

i=l,...,

/&jk-Si<)Ui

&k>O Si>O

i=l,...,m

Vi,j,k Vi

dj

Ilt, j=l,...,

n

(1)

(2)

(3) (4)

If we can find a dual solution that is feasible in the dual problem, and a primal solution that is feasible in the primal, such that the complementary slackness conditions are satisfied, we have found the optimal solution of the LP-relaxation of the location problem with spatial interaction. Furthermore, if the primal solution is integral, it is also the optimal solution of the simple plant location problem with spatial interaction. Now we note that if TVis fixed in the dual, it is easy to find a dual feasible solution and a primal solution that satisfies the complementary slackness conditions and also to determine whether the primal solution is feasible. If it is not, we also get indications of how to change a, in order to make it feasible. This forms an excellent basis for a dual search method in CC

Dual search procedures

91

Hence we will construct a dual ascent and adjustment method in a for solving uLp=max h(a) where h(a) is the dual function given in section 4. Here we note the similarities to dual subgradient optimization, as described in the previous section. Both are search methods for solving the dual, but will use very different techniques to search through the dual space of a.

6. THE DUAL ASCENT AND ADJUSTMENT METHOD To evaluate h(a) we assume that a is temporarily lixed to G. The dual subprobfem is then separable in i, in the following way.

h(ii)=f

V!(i)+i i=l

WjEj j=l

where, Vi, ~~(~)=ma~-d,

Vj,k and

s.t.-&<&k-&j

constraints 2, 3, 4 of the LP-dual.

In Holmberg and Jijmsten (1995a), we show that a dual feasible and optimal solution is obtained by Vi,j,k

&~=m~(O,&-E+) and

which yields

h(a)=i

Wjajj=l

f

max(O, Ctj-?ijk)

max j=l

i=l

k=l

. >

h(a) is clearly a concave, piecewise linear function. We will now maximize h(a) with a greedy heuristic, which is what Erlenkotter (1978) did for the ordinary uncapacitated facility location problem, and which gave the very efficient algorithm DUALOC. In the following, we consider a fixed to &. Definition: Let, for each (i,j), kijE{O,l,.. ., Wj} be 5UCh that &ZEijk Vkk, and let qije{O,1) be such that 1 if qij=

0

&j=Eijk,,

if not

K. HOLMBERG and K. J6RNSTEN

92

Also let, for convenience,

j=l k=l

j=l k=l

Then we cu define the index sets Z”={i:g>yai}, I=={i:Si=yai}, Z’={i:fi<. Holmberg and Jiimsten (1995a) show that the complementary slackness conditions can be reduced to the following. rijk=ri

Vk
(Cl)

x+=0

Vk>k,,Vi,Vj

((2

Vid’

(C3)

Zi=O ‘did’

(C4)

Zi=l

where the index sets can be easily calculated from any given Ccas shown above. Now we have a dual feasible solution that satisfies the complementary slackness conditions, so in order for this solution to be optimal, the primal solution should be feasible. Constraints 2, 3, 4 and 5 are trivially satisfied, by letting Zi attain any value between 0 and 1 and x&kany value between 0 and Zi. The difficulty now lies in satisfying constraint 1. Upper and lower bounds on the left-hand side of constraint 1 are obtained as follows, as shown in Holmberg and Jiirnsten (1995a). Lemma: Let, Vj, W;(E)= C

(k,-qi,)

and WY(E)= c

i&l’

kF

id>

Then

The primal solution is thus feasible if Vj

(11) So if inequalities 11 are satisfied, we have an optimal solution to the LP-relaxation. For an optimal Cr,we get certain w;(E) and WY(&)Vj, which are associated with two primal integer solutions, not necessarily feasible. If wi(&)=wy(5) Vj, we have one feasible, integer solution, and there is no duality gap. Our solution approach is to increase a in small steps, in order to satisfy the inequalities in Il. Suppose that we change one element in the a-vector, i.e. let a,=&!+&, where 0,>0, and aj=Gj Vj#I. We have Ei&,,I El<&,,+ 1 Vi, and we make the change small, so that Ci&,, C 5, + 0, ,
Dual search procedures

93

To identify the changes in the index sets Z, we note that s~=s~+~&. All elements in I= with kil>O will enter I’. By making the change small enough, we can make sure that no element goes directly from I’ to I’, but stays in I=. To achieve this, we calculate an upper bound on the change, denoted by 6,. If id’, i.e. Si< +tJUi, we require that Si+ kilOl
where

The bounds w:(G) and w;(G) Vji.are changed in the following way, for a certain change 6,. Let Z;l=(iEZ>, fil=Eilk,,+l -till and denote the change in w;I(cL)with AwY(tl) and the change in

w;(E) with Aw:(S). Then AwY(&)=lZ’iland Aw:(&)= 1

qib

is1’

AISO let Z!={id’:kg>O,

B,=k

(rai-fi)},

Z’;{ieZ’: kg>O, t$=-$

II gI+iIk,,

+

(y~i-._?i), rl

1-c(i) and Z;={id=:

kil>O}.

Then

kil+IZ’;I+lZt;I+

AwY(k)=C

ier;

1 kil+ 1 id;

qilVj#l

id’

and

AWN=

C id’;

k,+

1

(kij-qij)Vj

#I.

id;

Note that Awj(&)
K. HOLMBERG and K. JbRNSTEN

94

We only consider changing a certain aI as long as w?(i) < wI. When w;(i) awl, cclis optimally fired. We also refrain from increasing a certain aI if this makes wi(&)>wj for some j. In this case 4 is said to be blocked. We continue to increase &jin small steps, for each j, unless it is optimally fixed or blocked, until all variables are optimally fixed or blocked. Now let .ZF={j: w,!(a)wj for any i&F then let JP=JPU{j} and JF=JF- {I}. If WY<@’ for any jd, and if #$ Wj, let J,=J,-(j). 6. If JF # C$and J, # 4, go to 1. Otherwise: If JI=$: Stop. The optimum is found. If JF=$: Stop. All variables are blocked. Go on to the adjustment procedure. The selection of j in step 1 can be made in several ways. The most simple method is to choose cyclically over jd, n JF. Another possibility is to choose the j for which the distance from the feasible interval [wj,wy] to the demand is the largest, i.e. at each iteration choose j as 0. Calculate

In the dual adjustment procedure, we consider decreasing some al ~EJ, that blocks further improvement in the dual ascent procedure. This means that we temporarily allow WYto decrease, even if this increases the distance to wP This is done in order to allow increase of other aj’s so that the total effect is positive.

Dual search procedures

95

Jp contains the indices of the CQ’whose S decrease may improve the solution, We know that i.e. the blocking variables are optimally fixed. Since the increase of each a[ is blocked, we know that for each 1 there exists some id;, such that any increase of al would increase si and make i enter I’. There also exists some jEJp, j # 1, such that

J&J,,

W,!+

C

(kij-qij)>Wj

id;

We want to increase aI further. This can be enabled by decreasing ap In this way, it is always possible to get around a blocking situation. Assuming that the index sets I’, I= and I’ are unchanged, we decrease aj, step by step, until

C

(kij-qg)=wj

iel'

which is done by repeatedly changing aj to next lesser cqk, i.e. repeatedly decreasing kij by one in an a priori unknown sequence of i. We call the procedure of decreasing aj in this way DECREASEQ). After this, we return to the dual ascent procedure, but do not allow any increase of aj (until all the other @j’sare blocked). In other words, we enter the dual ascent procedure with JF= { 1,. . . , n}\O_}. Now we give the main framework for the whole procedure, including the dual ascent and the dual adjustment parts. 0. Let JF=JI={l ,..., n), Jp=O and J,=fL For Vj, let ij=E~l=m~EijlVj,

k,=l, q,=l and kij=O, qu=O Vi #t.

1. Do ASCENT. If J,=$, stop. 2. Chooseaj~~.LetJo=J,u{i),JI=JF={1,...,n}\JDandJ,=~. 3. Do DECREASE (j). 4. Do ASCENT. If JI=& stop. Otherwise go to 2. Steps 2-4 are repeated a predetermined number of times or until an optimal solution is found. This is called the DUAL ADJUSTMENT phase. When the optimal Ccis found, we calculate the rest of the optimal dual solution as shown in section 3. The primal solution is calculated from the complementary slackness conditions. After this, we can apply a primal heuristic of greedy type to find a good primal feasible integer solution.

7. COMPUTATIONAL RESULTS As test problems, we have used data from the two problems (size 5 x 8) solved by Erlenkotter (1978) and subsequently by Jacobsen (1988) and a number of larger random problems. The data is taken from Erlenkotter (1978) in the sense that the costs ci]Wjare taken from Erlenkotter, and we have extracted the largest even number from each column to get Wj. (One can note that the larger the Wj’s are, the larger and more difficult SPLPE becomes.) Prohibited assignments are indicated by cij=lOO. The random problems are generated according to the following model. Coordinates for the locations and the customers are exponentially distributed around two centers, and the transportation costs, I+, are given by the Euclidean distances between the points. The flxed

K. HOLMBERG and K. JORNSTEN

96

Table 1. Data for problems slex81 and slex82 (from Erlenkotter, 1978).c denotes the cost matrix, w the demand and a’ and a’ the two sets of fixed charges Problem data

C

w

12 21 18 21 17 10

18 100 19 19 1.5 10

10 15 11 1.5 11 10

100 16 13 12 10 15

12 11 10 13 14 5

100 14 100 8 13 15

18 11 100 16 20 10

100 11 13 8 100 15

a’

a2

100 70 60 110 80

200 200 200 400 300

costs, ai, and the demands, wj, are uniformly distributed, Wj in the interval 2-50, and ai in an interval depending on the total demand, n and m. One of the problems was duplicated and supplied with different sets of fixed costs. The parameter y, reflecting the relation between the transportation costs and the effects of the spatial interaction, is very important, both for ‘the structure of the solution and for the difficulty of the problem. The value of y must be chosen carefully, and verified by real life data. For large values of the y linear parts of the objective function dominate, and for very small values of the y nonlinear (logarithmic) parts of the objective function dominate. At this stage, it is not possible to determine y by calibrating the models against reality, so therefore we have, in the following comparisons, used several different values of y, in order to show the efficiency of the methods under different circumstances. However, in a real life situation, only one value of y is of interest. For the Erlenkotter examples (data given above), we have compared the solutions for different values of y. For very small values, y=O.OOOOl, the transportation costs hardly affect the solution, which means that the spatial interaction mechanism (the logarithmic part of the objective function) spreads out the trips (all costs are very low). All facilities are then opened. On the other hand, for large values, y=5.0, the logarithmic part of the objective function (the spatial interaction) is negligible, and the optimal solution is identical to that of SPLP. In Table 2, we give the problem sizes for the different problems (both m and n as well %s the number of variables, NV, and constraints, NC, of SPLPSE) and additional problem data. In Table 3 we give, for the different values of y, the CPU time (in seconds on a DEC system 5500) and the relative error E, obtained, for the dual ascent method and the subgradient method, on average for the different groups of problems. The programs are written in Fortran 77. More detailed computational results can be found in Holmberg and Jdrnsten (1994). We use three parameters to specify different strategies of the dual ascent method. The first one, pl, determines whether the dual variables mj should be increased cyclically (pl=O) or if the one with maximal residual should be chosen first (&= 1). The second parameter, p2, decides how much + should be decreased in the adjustment phase, and is set in one of three ways: (1) so that the lower bound coincides with the demand (p2=O), (2) one step more than in the previous case, i.e. so that the lower bound becomes less than the demand &=l), (3) or just a single step @=2). The third parameter,p3, controls whether (&=l) or not @=O) at most five additional attempts of the adjustment and ascent procedures should be done. In the tables of computational results, two parameter settings are included, DA1 with (0, 2, 0)

Dual search procedures

91

Table 2. Description of the test problems Problem

#P

m

n

N”

NC

KUAX

slex81 slex82 slex2 sleda slex2b

15 15 15 10 8 8 9 I 5 5 5 5 4 4 4

5 5 8 8 8 8 8 10 20 20 50 50 50 50 50

8 8 8 8 8 8 8 20 50 50 100 100 100 100 100

455 455 1464 1464 1464 1464 1464 4910 26300 26240 115OOb 115900 117000 116850 116050

458 458 1464 1464 1464 1464 1464 4920 26330 26210 115050 115950 117050 116900 116100

15 15 43 43 43 43 43 46 50 50 50 50 50 50 50

slex2d slex3 slex4 slexal slexa2 slexd

WlWT 90 90

182 182 182 182 182 490 1314 1311 2299 2317 2339 2336 2320

Comment Erlenkotter 1 Erlenkotter 2 Random ex slex2 with f= 10000 sled with f=5000 slex2 with f= 1000 slex2 with f= 100 Random ex Random ex Random ex Random ex Random ex Random ex Random ex Random ex

#P denotes the number of instances treated (i.e. the number of different values of y used), NV the number of variables, NC the number of constraints, wMMAx the maximal demand and wToT the total demand

and DA2 with (1, 0, 1). The subgradient method was run in at most 25 iterations (SG25) for the smaller problems and at most 10 iterations (SGlO) for the larger ones. The tables reveal that the dual ascent methods are superior to the subgradient method for small values of y, while for larger values of y, the subgradient method dominates. For the largest problems, the dual ascent method DA1 yields very poor bounds, while DA2 yields much better bounds, especially for small values of y. For large values of y, however, DA2 exhibits a significant increase in solution time. On the other hand, the subgradient method yields better bounds and even shorter solution times as y grows. The explanation for this is that the dual ascent method draws heavily on the logarithmic part of the objective function, by making the steps so small that no breakpoint of the dual function is bypassed. The sizes of these steps are mostly determined by the difference which is only due to the logarithmic part of the objective function. between cijk and For large values of y, these steps are insignificant, and thus it is not a good method. On the other hand, the subgradient method does not depend too much on the logarithmic part, and the step length is determined in a rather simplistic way. Looking for the “best” method, one can use different criteria. Considering only the accuracy, we find that DA1 yields the least accuracy in 34 cases, DA2 in 64 cases, SG25 in 88 cases and SGlO in 11 cases. (In some cases there are ties). A weakness of this criterion is shown by the fact that SGlO can only be best when SG25 is not used or SGlO and SG25 yield the same accuracy. Taking the solution time into account by looking for the method that gives a better accuracy in less time than any other method yields the following. In 21 cases DA1 is the best, in 34 cases DA2 is the best, in 6 cases SGlO is the best and in none of the cases SG25 is better than the other methods. A weakness of this criterion is indicated by the fact that, in many cases, no best method is found, since a better accuracy is found in a longer time (or a worse accuracy is found in a shorter time) by another method. All these numbers must be considered as very approximate, but we still draw the conclusion that DA2 seems to be more efficient than the other methods. The subgradient cijk+

1,

K. HOLMBERG

98

and K. JGRNSTEN

method is in most cases slower than the dual ascent method, but it has the advantage that a better solution can in many cases be found by simply allowing the method to continue for a larger number of iterations (say 100 or 200). The dual ascent method cannot be continued in such a simple way, this would require modifications of the adjustment phase. However, it is clear that the dual ascent method is quicker than the subgradient method for reaching the accuracy obtainable by both methods, which in many cases is satisfactory. In any case, it is interesting to note that the two methods complement each other very nicely. For small values of y, the dual ascent method is very efficient, while for large values of y the subgradient method is very efficient, and together they give possibilities of solving all of the test problems efficiently for virtually all values of y. We also note that the formulation SPLPSE seems to be a very strong one, in the sense that the duality gaps in general are very small. This has been verified by the solution of the smaller instances (with not more than 5000 variables) by an Integer Programming code. In many of these cases, an integer solution was immediately found; the difference between uLp and ZI* was never larger than 1.5% and the branch-and-bound tree never contained more than 12 nodes. Table 3. Average solution times for the dual ascent method (two parameter settings) and the dual subgradient method (at most 2.5 or 10 iterations) Problem data

DA1 (0,2,0)

DA2 (l,O,l)

Name

y

Err

Time

Err

slex8*

o.txlOO1 0.0001 0.001 0.01 0.05 0.1 0.5 1.0 2.0 5.0

0.00 0.00 0.00 0.00 0.01 0.02 0.06 0.06 0.02 0.03

0.02 0.02 0.02 0.02 0.03 0.04 :: 0:os 0.06

0.00 0.00 0.00 0.00 0.00 0.02 0.09 0.04 0.12 0.11

0.00001 0.0001 0.001 0.01 0.1 1.0

0.00 0.00 0.06 0.05 0.08 0.06

0.07 0.08 0.08 0.18 0.26 0.25

0.00001 0.0001 0.001 0.01 0.1

0.02 0.87 0.22 0.12 0.04

slex4 *

0.00001 0.0001 0.001 0.01 0.1

slexa*

0.00001 0.0901 0.001 0.01

slex2’

slex3

Time 0.03 0.02 0.02 0.02 0.05 0.10

Subgrad Err 0.00 0.00 0.00 0.00 0.00 0.01 0.07

SubgradlO

Time 0.07 0.08 0.12 0.11 0.08 0.10

Err

Time

-

-

:: 0:08 0.07

83 0:OO

:: 0:10 0.12

:Fz 0:06 0.06 0.06

0.06 0.06 0.08 0.98 1.69 0.79

0.00 0.00 0.00 0.00 0.10 0.10

0.46 0.41 0.26 0.28 0.27 0.27

-

0.26 0.14 0.20 0.35 0.63

:: Ok0 0.00 0.00

0.31 0.19 0.24 1.66 2.52

:: 0:OO 0.00 0.00

1.68 1.43 1.42 0.99 0.78

0.01 0.01 0.00 0.01 0.01

0.59 0.59 0.55 0.41 0.37

0.50 0.92 0.56 0.18 0.06

1.10 0.88 0.90 1.67 3.07

0.00 0.00 0.21 0.02 0.06

1.10 0.94 1.33 7.56 10.62

0.00 0.00 0.00 0.00 0.02

6.56 6.78 5.88 4.16 3.80

0.02 0.01 0.00 0.00 0.05

2.74 2.72 2.66 2.02 1.57

0.41 0.97 0.84 0.30

4.69

0.00 0.01 0.21 0.14

3.80 3.97 6.27 29.17

0.09 0.04 0.00 0.00

22.28 22.23 21.72 18.18

0.16 0.13 0.00 0.02

9.14 7.28 8.94 7.84

0.00

Dual, search procedures

99

An attempt to compare our solution methods to other solution methods found in the literature yields the following. Jacobsen (1988) solved problems of sizes 5 x 8 and 8 x 8 with a branch-and-bound method. The larger of these problems was solved with 12 different sets of fixed costs. Erlenkotter and Leonardi (1985) solved problems of sizes 8 x 8, 23 x 23 and 44 x 69. No solution times are given in either of these reports, but comments by Erlenkotter and Leonardi indicate that their solution times are not very short, and that problems with large fixed costs are difficult to solve. In their conclusions, they note that dual based solution methods are very promising, but that no natural extension of Erlenkotter’s method for location problems with spatial interaction seems evident. The dual ascent method presented in our paper is an extension of Erlenkotter’s method (though maybe not so evident), capable of solving larger problems than were previously solved. The limitations of this method, occurring for larger values of y, are overcome by combining it with the subgradient optimization method. 8. CONCLUSIONS We present a new “exact” formulation of the simple plant location problem with spatial interaction, which combines the costs for transportation and opening the plants with a spatial interaction term in the objective function, without using any approximation. An exact linearization of the somewhat awkward objective function yields a linear O/l-programming problem, with a structure similar to the ordinary simple plant location problem. Each transportation cost coefficient in the model consists of a logarithmic part (from the spatial interaction term) weighed together with the original transportation cost coefficient. We believe that the model can be very useful in practice, in situations when the public can choose the facility they patronize, and we also show that the formulation has advantages when it comes to solution methods. We present two dual search heuristic procedures for solving the problem. One is based on dual ascent and the other on subgradient optimization. As the computational tests show, the methods are very efficient when used for different values of y. Using the dual ascent method if the value of y is small and the subgradient method if the value of y is large produces a quite efficient solution procedure, capable of solving large problems (with around 116,000 variables) to an accuracy of better than 1% in a quite reasonable time (a matter of seconds). REFERENCES Bazaraa, M. S. & Goode, J. J. (1979) A survey of various tactics for generating Lagrangean multipliers in the context of Lagrangean duality. European Journal of Operational Research, 32,322-338. Beaumont, J. R. (1980) Spatial interaction models and the location-allocation problem. Journal of Regional Science, 20,37-50. Bilde, 0. & Krarup, J. (1977) Sharp lower bounds and efficient algorithms for the simple plant location problem. Annals of Discrete Mathematics, 1,79-97. Coelho, J. D. (1980) A locational surplus maximization approach to public facility location. Nota 15, Centro de Estasica e Aplicacoes, Dept. Matematica Aplicada Faculdade de Ciences de Lisboa, Lisbon, Portugal. Coelho, J. D. & Wilson, A. G. (1976) The optimum location and size of shopping centres. Regional Studies, 4, 413-421. Erlenkotter, D. (1978) A dual-based procedure for uncapacitated facility location. OperationsResearch, 26, 992-1009. Erlenkotter, D. & Leonardi, G. (1985) Facility location with spatially interactive travel behaviour. Sistemi Urbani, 1,29-41. Fisher, M. L. (1981) The Lagrangean relaxation method for solving integer programming problems. Management Science,21, l-18.

100

K. HOLMBERG and K. JGRNSTEN

Geoffrion, A. M. (1974) Lagrangean relaxation for integer programming. Mathematical Programming Study, 2, 82-114. Geoffrion, A. M. & McBride, R. (1978) Lagrangean relaxation applied to capacitated facility location problems. AIEE Transactions, l&40-47.

Held, M., Wolfe, P. & Crowder, H. (1974) Validation of subgradient optimization. Mathematical Progrumming, 6, 62-68.

Holmberg, K. & Jornsten, K. 0. (1989) Exact methods for gravity trip-distribution

models. Environment and

Planning A, 21, 81-97.

Holmberg, K. & Jdrnsten, K. 0. (1990) An exact formulation of the simple plant location problem with spatial interaction. Proceedings of Meeting V of the EURO Working Group on Locational Analysis, Orban-Ferauge, Fr. & Rasson, J. P. (Eds). Namur, Faculte Universitaries Notre-Dame de-la-Paix, Belgium. Holmberg, K. & Jdrnsten, K. 0. (1991) The exact formulation of the simple plant location problem with spatial interaction and its solution. Working Paper LiTH-MAT/OPT-WP-91-09, Linkiiping Institute of Technology, Sweden. Holmberg, K. & Jornsten, K 0. (1994) Dual search procedures for the exact formulation of the simple plant location problem with spatial interaction. Research Report LiTH-MAT-R-94-38, Linkiiping Institute of Technology, Sweden. Holmberg, K. & Jornsten, K. 0. (1995a) A dual ascent procedure for the exact formulation of the simple plant location problem with spatial interaction. Optimization, 36, 139-152. Holmberg, K. t Jornsten, K. 0. (1995b) Decomposition methods for the exact formulation of the simple plant location problem with spatial interaction. Research Report LiTH-MAT-R-95-13, Linkdping Institute of Technology, Sweden. _ Jacobsen, S. K. (1986) Plant location by entropy maximization. Research Report 10/86, IMSOR-DTH, Lyngby, Denmark. Jacobsen, S. K. (1988) On the solution of an entropy maximizing location model by submodularity. Research Report 6/88, IMSOR, DTH, Lyngby, Denmark. Kiirkel, M. (1989) On the exact solution of large-scale simple plant location problems. European Journal of Operational Research, 39, 157-173.

Leonardi, G. (1978) Optimum facility location by accessibility maximizing. Environment and Planning A, 11, 1287-1305. Leonardi, G. (1983) The use of random-utility theory in building location-allocation models. In Locution Analysis of Public Building Facilities. J. F. Thisse & H. G. Zoller (eds). Amsterdam, North-Holland Publishing Company, 357-383. Poljak, B. T. (1967) A general method of solving extremum problems. Soviet Mathematics Doklady, 8,593-597. Poljak, B. T. (1969) Minimization of unsmooth functionals. USSR Computational Mathematics and Mathematical Physics, 9, 14-29.

Shapiro, J. F. (1979) A survey of Lagrangean techniques for discrete optimization. Annals of Discrete Mathematics, 5, 113-138. Shor, N. Z. (1985) Minimization Methods for Nondifferentiable Functions. Berlin, Springer-Verlag. Wilson, A. G. (1967) A statistical theory of spatial distribution models. Tmnsportation Research, 1, 253-269.