The obnoxious p facility network location problem with facility interaction

EUROPEAN JOURNAL OF OPERATIONAL RESEARCH ELSEVIER

European Journal of Operational Research 102 (1997) 302-319

The obnoxious p facility network location problem with facility interaction S.B. Welch, S. Salhi* Management Maths Group, School of Mathematics and Statistics, University of Birmingham, Birmingham, UK

Abstract Three heuristics are proposed to solve the maximin formulation for siting p facilities on a network considering a pollution dispersion equation and facility interaction. Initially, the single facility problem is approached by building up polygons which model pollution spread about the nodes of the network. This is extended in the first heuristic to the p facility problem. The second combines both the p-maximin and p-maxisum objectives in a lexicographic manner. The third is based on recent developments of Simulated Annealing. The proposed heuristics are evaluated for up to six facilities on a set of randomly generated networks having 20 to 40 nodes and low or medium arc density. (~) 1997 Elsevier Science B.V. Keywords: Location; Obnoxious facility; Heuristics; Simulated annealing

1. Introduction The traditional location problem is concerned with the location of one or more facilities, in some solution space, so as to optimise some prespecified criteria. Such facilities can be classified as: desirable, semi-obnoxious or obnoxious (undesirable). Most services are provided by desirable or non-obnoxious facilities. These may be supermarkets, warehouses, shops, garages or banks for example, where it is beneficial for the facilities to be sited close to the customers that they will be serving. Sometimes, a facility produces a negative or undesirable effect, which may be present even though a high degree o f accessibility is required to the facility. If the undesirable effect outweighs the accessibility required, then the facility is said to be obnoxious. Some examples are: nuclear power stations, military installations and pollution producing industrial plants. Although necessary for society, these facilities are undesirable and often dangerous to their surroundings. Urban areas may endure lower local house prices and quality of life, whilst areas of natural beauty, e.g. National Parks, rivers and forests, may be damaged irrevocably. To minimise these effects the location of obnoxious facilities within an environment needs to be carefully addressed. In the non-obnoxious set of location problems, the objective function to be optimised is often of a pminisum or p-minimax form. For the semi-obnoxious problem the above objectives are often combined with a * Corresponding author. 0377-2217/97/$17.00 (~) 1997 Elsevier Science B.V. All rights reserved. PII S 0 3 7 7 - 2 2 1 7 ( 9 7 ) 0 0 1 1 1 - 2

S.B. Welch, S. Salhi/ European Journal of Operational Research 102 (1997) 302-319

303

lower bound on the distance between population centre and facility; see Moon and Chaudhry [18], Moon and Papayanopoulos [19] and Fonseca and Captivo [7]. In the case of obnoxious facilities the problem is reformulated as: (a) p - m a x i a n / p - m a x i s u m problem,

(b) p-maximin problem.

The p-maxisum problem (a) is linked to the p-minisum problem, in that the sum of the distances is maximised rather than minimised. Similarly, the p-maxian problem is related to the p-median problem in that its objective is to maximise the sum of the weighted distances between population centres and their nearest obnoxious facility. Church and Garfinkel [2] solved the l-maxisum problem on a network using shortest path distances in O(m2), where m is the number of nodes. Ting [25] improved this to O ( m ) for tree networks. In (b), Drezner and Wesolowsky [5] solved the single facility planar location problem by drawing a circle around each point and taking the area outside their union whilst satisfying some upper bound constraints. A type of bisection method was then used to locate the optimal solution. This concept will be extended in this study. Dasarathy and White [ 3] suggested a method based on the Kuhn-Tucker conditions for the unweighted problem of locating a point in a given convex polyhedron. Melachrinoudis and Cullinane (MC) [ 16] also used the Kuhn-Tucker conditions for the maximin problem, where each node is assumed to occupy a region which is approximated by a circle. Karkazis and Karagiorgis [ 13] extended the previous formulation to a region containing a set of areas defined as either protected or forbidden and not necessarily circular in nature. The problem of siting a facility becomes more difficult when considering pollution dispersion subject to meteorological conditions present. Karkazis and Papadimitriou [ 14] considered this problem in the plane. The objective proposed was the minimisation of the sum of the weighted risk factors for each node summed over all possible wind directions. A similar technique for the network was presented by Karkazis and Boffey [ 10] and Karkazis, Boffey and Malevris [11 ] who recommend using an interpolation scheme to estimate the data between the discrete wind directions. The aim of this study is to: (i) develop a heuristic to site a single facility on a network satisfying the maximin criteria, which approximates the pollution dispersion by the use of polygons; (ii) investigate the p facility problem, where the facilities interact in some way, subject to pollution dispersion; (iii) combine both the p-maximin and p-maxisum objectives using a lexicographic approach; (iv) produce Simulated Annealing implementations to assess the performance of the proposed heuristics. We describe (i) in Section 2, (ii) in Section 3.1, (iii) in Section 3.2, (iv) in Section 3.3 and present our computational results in Section 4. We close the paper with conclusions and some suggestions.

2. The single facility maximin problem A heuristic procedure is first introduced for solving the single facility network location problem. This will serve as a basis for the development of the heuristics which locate p facilities discussed in the next section. The network involved denotes a transport system which may be road, rail, or a set of flight or shipping routes. The nodes represent a set of existing bodies within the network that may affect the location of the facility. These may be population centres, existing facilities or environmental bodies. It is assumed that the facility and nodes of the network can each be approximated by a point and the arcs of the network are straight line links. Formally, the problem is to site a single facility X on a network G = ( N , A ) , consisting of a set N = {Xl . . . . . Nm} of nodes and a set of arcs ( i , j ) C a with i , j E {1 . . . . . m} so as to maximise the minimum weighted Euclidean distance from the facility to the set of nodes. This can be formulated as

S.B. Welch, S. Salhi/European Journal of Operational Research 102 (1997)302-319

304

max I f ( X ) = minwid(X, Ni)2, / = 1 . . . . . m], XCG L

i

(1)

J

where wi is a non-negative weighting of node Ni and d(X, N/)2 is the Euclidean distance between facility X and node Ni. In this section, we propose an algorithm which solves Eq. (1) by finding the feasible location intervals within a prescribed error bound. For any given weighted distance wd, a series of circles ci of radius wd/wi can be formed around the existing nodes Ni of the network, where each point on a circle is of equi-pollution to all others. The points of equal or higher pollution are defined by the convex hull of ci, conv(ci). The arc lengths FA containing those arc points that do not lie within any convex hull (FA = G - Ui conv(ci)), correspond to the feasible solution intervals for the given weighted distance. An iterative procedure is then evoked that converges towards a weighted distance lying within a prespecified bound. The formulation presented above contains one major underlying assumption that may restrict its applicability. It assumes that the amount of pollution at a point is inversely proportional to the square of the distance from the facility to that point. That is, the spread of pollution is circular in nature and independent of all environmental occurrences. This is very limiting and may provide a solution that is misleading. To overcome this, we proceed to a more accurate pollution dispersion model by considering a dispersion equation, allowing a more intimate study of the relationship between pollution levels, distance and wind strength. One such frequently used dispersion equation is given by Karkazis and Papadimitriou [ 14] (see Appendix A for a detailed formulation). Using this equation, dispersion maps Di each homeomorphic to a circle replace the circles used previously. The problem therefore becomes finding those intervals outside the union of all maps present. As all points within a map Di are point-convex with respect to node Ni, then the arc lengths FA can be defined as FA = G - UPconvN,(Di), i

where

p c o n v N , ( D i ) = ~AX + (1 -- /~)Ni : 0 <~ A <~ 1, x C D i ~ . )

By studying the dispersion equation, over the active interval, it can be shown that contrary to the previous assumption the pollution concentration monotonically increases up to a point rh and then monotonically decreases. This causes problems to occur in the solving by the method previously established. Each dispersion map represented a contour line of equi-pollution concentration, where all points within, represented by pconvN, (Di) - Di, have a greater concentration and similarly all points outside are at a lower concentration. This is no longer the case and in order to solve this problem two maps must be formed, say D a and D b. Both represent equi-pollution concentration contours, but only those points that lie within D = pconvui(Da) -- pconvui(D b) -- D a are at a higher pollution concentration. The problem is therefore reduced to solving two sets of dispersion maps for each node such that the infeasible region is made up of the set of points between each pair of maps. However, this raises the issue of applicability. The above two map method may locate the facility next to a node of small weighting (small implies important in this context). This may be ideal in the continuous emission case but when other factors are considered, such as accidental release or public opposition, could be disastrous. To prevent this 'close' siting, a lower bound should be imposed around each node of the network. This bounding may be implemented in the form of a straight inequality being satisfied or by the use of a maxi/minisum formulation as a secondary objective. For the purposes of the algorithm to be presented a lower bound of distance rh is imposed. This simplifies the problem, but for large values of r h may be impractical. In such circumstances it may be necessary to instigate the use of the two map method, with lower bounding. The formulation obtained is similar to that previously stated, in that max

XEG

[f(r,O)=

miinwi/Ci(r,h,O,O), / = 1 . . . . . m]

s.t. d(X, Ni)2 >>-rh

Vi,

(2)

S.B. Welch. S. Salhi/European Journal of Operational Research 102 (1997) 302-319

305

where Ci(r, h, 0, 0) is the pollution concentration at distance r and direction 0. As stated in Karkazis, Boffey and Malevris [ 1 1] the data available for the building up of any wind speed rose using wind velocities and probabilities is often limited to only eight directions (N, NE, E, SE, S, SW, W, NW). This causes a problem in that the data is required to be continuous. To overcome this a clamped cubic spline interpolation technique is used between the discretely defined wind directions. 2.1. Basic idea of the heuristic

Combining the wind speed cubics with the dispersion equation implies that finding the relevant intersections with the arcs of the network would involve using some root finding method. This would take a large amount of computing time. Therefore, a polygon approximation of the wind rose is used. The error involved in this approximation is defined by both a perpendicular error between the sides of the polygon and the cubic, and a radial error from the node about which the polygon is centred. This error is decreased progressively through the algorithm by increasing the number of vertices defining the polygons. Now for any given weighted inverse pollution level wp of the form given in Eq. (2), the wind velocity and pollution rate of emission are inputed into the dispersion equation which is then solved to give a radial distance r for the given direction. This distance must be calculated by a root finding method. Several methods were tested, but care needs to be taken to prevent any given iteration producing a solution outside the defined interval. This may cause the exponential part of the equation to be zero, causing the iterative process to break down. The error is then calculated between each side of the polygon P and its corresponding interpolated dispersion polynomial D and transformed into a dispersion error. Since those intervals, FA, outside [-JiPc°nvo,(Pi) are feasible, subject to the polygon approximation error, then an iterative procedure is evoked to find a feasible facility location within an initially defined error bound. The formal algorithm, Dispfloc, is given in Appendix C. Within the algorithm special attention should be paid to deltastep. This is the size of the change in the weighted inverse pollution level, wp, at any given iteration. Careful monitoring of the error and the change in deltastep is required to avoid the creation of a lower bound for wp that is infeasible and cannot be crossed due to the continual decrease in the size of deltastep - see Welch and Salhi [24] for more information.

3. The p facility problem with facility interaction Often, several facilities must be sited where the facilities interact in some way. In the obnoxious case it may be preferential to site all the facilities as far away from each other as possible. For example, in the case of military installations to minimise the attack potential of opposing forces. This section introduces the p facility location problem, where it is necessary for the weighted facility interaction to also satisfy the maximin criteria. A constructive heuristic based on that given previously is proposed. The result of this p-maximin formulation is then set as a constraint to a p-maxisum approach, which is again solved by a constructive heuristic. The results of this are then compared with those obtained by three variants of Simulated Annealing. 3.1. The p-maximin location problem

This idea of facility interaction causes the formulation of the problem encountered in the previous section to alter slightly. The problem is to site p facilities XJ, J = 1. . . . . p, on a network G = ( N , A ) , consisting of a set N = {N~ . . . . . N,,} of nodes and a set of arcs ( i , j ) C A with i , j E {1 . . . . . m} so as to minimise the maximum weighted pollution dispersion level. This is then transformed to that of a maximin form to give the following formulation:

S.B. Welch, S. Salhi/European Journal of Operational Research 102 (1997) 302-319

306 max

XjEG s.t.

[G(r,O) = m i n { i.j

wi Cx.u,(r,h,O,O)

uj }] ' Cx, xj(r,h,O,O) ,

i= l . . . . . m,

Cx, ui(r,h,O,O) = m a x {Cx,,N~(r,h,O,O) . . . . . Cxj,,u,(r,h,O,O)},

j= l ..... p

i= 1 . . . . . m,

(3) (4)

Cx, x~ (r, h, 0, 0) = max {Cx,,xj (r, h, O, O) . . . . . Cxj ~,xj (r, h, O, 0), Cxj~,,xj(r,h,O,O) . . . . . Cx,,,xj(r,h,O,O) },

j = 1 . . . . . p,

(5)

where Cx,,xj (r, h,0, 0) is the pollution effect from facility i acting on facility j, vj is the facility interaction weighting and all other notations are as previously. To take into account the new facility interaction, several changes must be made to the algorithm Dispfloc. The error boundings used so far have been upper bounds on the radial wind speed error, ew, and radial distance error, ed, respectively. No distinction has been made, in either of these two boundings, between a positive and negative error. In the p facility case, this lack of segregation has important connotations. At any given instant within the algorithm Dispfloc, the radius of the given interpolated wind pollution level can be below or above the polygon radius by ed. This implies that an infeasible solution could be made feasible for the p facility case by increasing the number of vertices defining the wind velocity polygon and so decreasing the error. This was overcome in the single facility case by the careful monitoring of both the error bound and the feasible/infeasible arcs. In the p facility case, this monitoring will be extremely difficult. Not only this, but any small miscalculation in this monitoring could lead to large changes in the p facility solution. It seems sensible therefore, to separate the error-below, edb, from the error-above, ed.~. A list of arcs containing a feasible interval irrespective of the interpolation accuracy are obtained. Also, a set of feasible solution intervals FA = f(17•(i'J) l~.~Lk, 1 , lr~A(i'J) . . k , 2 ), ( i , j ) C A, k = 1 ' ' " .} are given that do alter with the interpolation accuracy. This enables the algorithm to distinguish between an infeasible solution that will always be infeasible for the given wp and a solution that would become feasible by increasing the accuracy of the interpolated pollution curve. The stopping criterion used previously is now invalidated by the fact that within any feasible solution for the p facility problem, the maximum length of any interval in FA does not lie within a specified tolerance range. To replace this, a method that bounds both deltastep and the maximum error em (linking ej,,e~ h, the node weightings and the angle subtended by the node and feasible intervals) is used, so that the solution is obtained when both deltastep and em are sufficiently small. This sounds adequate but may result in the algorithm stopping prematurely. If the error bound has already been satisfied within the computation, then by the above stopping criteria all that is necessary is to reduce deltastep so that it also satisfies the above conditions. This will result in an inaccurate solution. To prevent this, a third stopping criterion is introduced that states that the previous iteration must have given an infeasible solution. Essentially the p facility interaction problem is solvable in exactly the same manner as the single facility problem, the difficulty being that the sites at which the polygons are centred (each of the p facilities, as well as the nodes present) are unknown. After choosing the p sites for the facilities the algorithm can be resumed and the feasibility of the iteration in question confirmed or denied. The difficulty becomes apparent when it is realised that the question posed is not: "Is one particular siting infeasible or not for wp?", but is: "Are all possible sitings of the p facilities infeasible for any given wp?". Here a step-by-step placement of the facilities is applied. Each facility will be sited in turn and the infeasible parts of the so far feasible solution intervals FA removed. This will continue until all facilities have been sited or it has been proven that for the given value of wp there exists no feasible siting solution. Two questions are then raised. First, in what order are the facilities to be sited and secondly, where in the continuous feasible set are they to be sited? Both of these questions are solved by introducing a measure/function of the 'sitability' of any of the facilities

S.B. Welch, S. Salhi/European Journal of Operational Research 102 (1997) 302-319

307

at any given point on the network. This function is based on the idea that feasibility has a greater probability of being attained if, after each iteration, the largest possible length of feasible solution intervals FA remains. In contradiction to this though, the facilities should be sited in order of increasing weight (decreasing importance), as the lower weighted facilities are more likely to affect the solution set and thus feasibility. This ordering of the facility placement may not be a definitive approach to this problem, as in certain cases the results obtained may be misleading. If this ordering is at any time deemed unsatisfactory it can easily be removed by solving the problem for each permutation of the step-by-step process. The step-by-step placement may also cause discrepancies in the results compared to a completely random siting of all p facilities at one time. However, it does have the advantage of being quicker than the above and far more guidable. This step-by-step method would be fairly straight forward in the discrete case, but unfortunately here each facility can be sited anywhere on the continuous feasible intervals FA, found at each iteration, thus complicating the problem. This difficulty is overcome by assuming that if a facility can be located along a feasible interval FA, then it can be sited feasibly at either of its end points. This is preferential to an 'interior' location by considering the previous statement of minimising the feasible interval lengths that become infeasible at each step of the method. If one considers only the interval in question, then it is obvious that end placement is beneficial. If one includes all other intervals present, it no longer remains clear. To aid in this, the function below is introduced, as a measure of the 'closeness' of intervals to the end point in question. (a)

r

(b)^

1

(i j) f(FAk, / ) = A;FA(i,i)FA(~,.I~ ~

"~k,1 '

k,2

(,,,,)6a ~kw=, ~-'~2z=,d(FA ~i),FA~Uz*,))2 )2q-a

+

/~*P

,

(6)

where (i,j) E A; I = 1,2; k = 1 . . . . . no. of solutions on ( i , j ) ; and a, /3 are parameters. This function is dependent not only on the length of the feasible interval upon which the endpoint is to be sited (part ( a ) ) , but the sum of the distances to all the other endpoints of the feasible intervals FA ( p a r t ( b ) ) . Formally the main steps of our algorithm, which we refer to as Facloc, are given in Appendix C. The algorithm Facloc can now be combined with the previously stated single facility location algorithm, Dispfloc, to give an algorithm that sites p facilities on a network where the facilities interact in some way. The main steps of this algorithm which is denoted by P-dfloc are given in Appendix C.

3.2. A lexicographic approach for the bi-objective problem It can be observed in Facloc that often a number of the facilities can be relocated elsewhere in the set FA, without violating the feasibility constraint or reducing wp. This implies that certain facilities may be resited to improve the solution of a p-maxisum formulation, whilst maintaining the bound set by the p-maximin problem. The new aim is to site p facilities on a given network satisfying the formulation of Eqs. ( 3 ) - ( 5 ) . When this is achieved, the facility that is causing the solutions to become infeasible for larger values of wp is sited and fixed. The algorithm then continues for the remaining p - 1 facilities until wp satisfies Eqs. ( 3 ) - ( 5 ) for p - 1 facilities. This is repeated until no further optimisation can be performed. This tries to overcome the difficulty of the maximin formulation, which only guarantees that all facilities are outside the minimum distance stipulated (in this case wp) and provides no measure of how far each facility is outside this bounding. This formulation strays more towards the p-maxisum idea, with the added advantage of prescribing a lower bound for wp. The problem is to site p facilities Xi, J = 1 . . . . . p, on a network G = (N,A), consisting of a set of arcs (i, j) C A with i,j E { 1 . . . . . m ) so as to maximise the minimum weighted inverse pollution dispersion level. This is repeated for p - k facilities for k = 0 . . . . . p - 1 until the optimal siting of all p facilities has been achieved.

S.B. Welch, S. Salhi/European Journal of Operational Research 102 (1997) 302-319

308

gp-k(r,O) = max [Gp_k(r,O)=min{ XjGG

i,j

wi vj }] Cx.ue(r,h,O,O)' Cx, xj(r,h,O,O) i = 1 . . . . . m,

s.t.

gn(r,O) = m a x Gn(r,O),

n=p-k+l

'

j = 1. . . . . p,

k=0 ..... p-I

. . . . . p,

Cx, Ni(r,h,O,O) = m a x {Cx,,N,(r,h,O,O) . . . . . Cx,,,N~(r,h,O,O)},

(7) (8)

i= 1. . . . . m,

(9)

Cx, xj ( r, h, 0, 0) = max { Cx~,xj ( r, h, O, O) . . . . . Cxj_, ,xi ( r, h, O, 0), Cxj~,,xj(r,h,O,O) . . . . . Cx,,,xj(r,h,O,O)},

j = 1. . . . . p,

(10)

where the same notation is used as previously. The formulation stated above can be solved using the existing algorithm P-dfloc. The only changes that need to be instigated are linked with the fixing of the (k ÷ 1)th facility when p - k facilities have been optimised and the effect those already fixed have on the feasible set of solution arcs FA. In fact these alterations can be introduced into P-dfloc by the introduction of several substeps. These can be observed in the algorithm P-subdfloc, given in Appendix C. An illustration of the heuristic P-subdfloc can be seen in Figs. la-d. These correspond to a network with m = 23 (nodes), p = 0.2 (arc density) and p = 3. Special attention should be paid to the change in the polygon approximation between Fig. la and lb and the fixing of a facility, reducing the problem to two facilities in lc. The solution is displayed in Fig. ld.

3.3. A Simulated Annealing approach Up until now, it has been assumed that each feasible interval FA can be approximated by a finite set of points. In the algorithms, P-dfloc and P-subdfloc, this took the form of limiting the set to only the end points of the arcs FA. Rather than using this discretisation, a Simulated Annealing technique can be employed, which retains the entire length of each interval FA. Simulated Annealing (SA) is a hill climbing heuristic which has the flexibility of accepting non-improving moves with a certain probability. This is contrary to local search descent methods (constructive heuristics) which accept only moves with positive improvement. SA has its origins in statistical mechanics, and Kirkpatrick et al. [15] proposed an SA algorithm based on the analogy between the annealing process of cooling solids and optimisation problems. For an in-depth guide to Simulated Annealing and other meta-heuristics (also known as modern or global heuristics) the reader is referred to Reeves [22]. The Simulated Annealing process will be used to investigate the efficiency of the results gained by the algorithm P-subdfloc. Three implementations of our Simulated Annealing heuristic will be used with which to compare the previous heuristics results. * The values for the cooling schedules corresponding to the algorithm SA1 are discussed below. No. of iterations: M = 6000. This provides enough iterations to allow the algorithm to escape from any large valleys of local optima, but not too many as to cause the computing time to be very large. Starting point: The result obtained from P-dfloc is used as a starting value. This is so the final solution is independent of the result gained by P-subdfloc. Initial temperature: To = --Amax/lnp, where p = 0.8 and Amax is the maximum {IF(s') - F ( s " ) I with s' and s H C N ( x ) } , the neighbourhood of x. Final temperature: Tf = 0.001. Obtained by observation. Neighbourhood: At any iteration the neighbourhood is variable and can be either of a shift type: where a facility is randomly chosen and resited anywhere on the feasible arcs that remain after its removal, or of an 4 exchange type: where two randomly chosen facilities swap location. The first is chosen with probability .~

S.B. Welch, S. Salhi/European Journal of Operational Research 102 (1997) 302-319

(a) : I t e r a t i o n

(c) : I t e r a t i o n

2

(b) : Iteration

15

(d) :Iteration

309

7

1 18

KEY • nodes facilities :::~facility interaction

• feasible interval : infeasible interval ::~ violated a r e a

Fig. 1. An illustration of P-subdfloc. and the second with probability ½. Temperature increment: The temperature at iteration k + 1 is updated using Eqs. (l l) and (12) depending on whether a move is accepted or not respectively; see Dowsland [4].

h+~ = Tk/ (1 + y r ~ ) ,

( il)

Tk+l = Tk/( 1 - o~Tk),

(12)

where/3 = A / M , ce = / 3 / 1 0 0 . 0 and A = (To - T f ) / ( T o . Tf). Stopping criteria: (i) Set at no improvement for 420 iterations; (ii) reaching the maximum number of

310

S.B. Welch,S. Salhi/European Journal of OperationalResearch 102 (1997) 302-319

iterations set at 6000; (iii) violating the lower bound on the temperature, Tf, whichever condition is violated first. * The second set of results, denoted by SA2, are obtained using the parameters defined below: Change in temperature increment: The temperature is decreased when accepting a move, according to Eq. ( 11 ), but the definition of/3 is now redefined as A ilk-M/o~+yx/~'

where y = 0 . 5

and o~= 10.0.

(13)

In contrast to previously, the temperature is only increased after 80 consecutive rejections of uphill moves. When this occurs, the temperature is reset to min{7~, 0 . 4 , TR} (TR is the previous reset temperature, initially TR = To); see Osman and Christofides [20] for a similar resetting scheme. Stopping criteria: A maximum number of 5 resets is performed without improvement, the final temperature is achieved or the maximum number of iterations set at 6000 is reached, whichever occurs first. • Finally, a post-optimiser, SA3, is used where all parameters chosen are as for SA2 but the starting point is the result obtained by P-subdfloc and not P-dfloc, as before. The Simulated Annealing process is used here: (i) to enable a direct comparison of performance between SA1, SA2 and P-subdfloc; (ii) to obtain best upper-bounds (SA3) with which to compare the performance of P-subdfloc, SA1 and SA2.

4. Computational results Our proposed heuristics are written in Fortran 77, using Uniras 6v4a for the graphical output and run on a Sun SPARC 20 workstation (712MP dual processor, 70/75 MHz with integer SPECrate 5726 and floating point SPECrate 5439). The algorithms described are tested on networks having 20 to 120 nodes (m) and arc density (p) of 0.2, 0.5 and 0.8 for the single facility case and for p facilities m = 20 to 40 and p = 0.2 or 0.5. Within each category, 20 networks are tested using Q = 100 g/s and h = 20 m. Our results will be based on the mean value of wp over these networks. We generated our networks as follows: the weights wi and the co-ordinates (xi, Yi) for each node Ni are generated randomly over {0 ..... 10} and (0, 125) respectively. The number of arcs incident to each node Ni is also randomly generated between ( p - 0.2) • ( m - 1) and p * ( m - 1), such that 80% of the arcs are within 0.4 • (maximum distance) to any other node and 100% of the arcs within 0.6 * (maximum distance). If the total number of arcs incident to each node cannot be satisfied, then the remaining arcs are ignored. In all of the networks tested, the solution accuracy is set to 200 metres. Given that the facility is sited half way along this feasible segment, then when the positive and negative polygon errors are taken into account, the facility is within a worst deviation of 200 metres. We first produce the results of the single facility problem obtained by Dispfloc in terms of computing time. The solutions of this algorithm are not reported, as these fall within an accuracy of 200 m of the optimal solutions. However, the computing times may be used as references for further comparison. The results found for the p facility cases are analysed using both solution quality and computing effort.

4.1. Single facility results Several parameters were highlighted within Dispfloc, upon which the efficiency of the algorithm depends. These were obtained by empirical testing: Yl = 2.0; Y2 = 2.0; Y3 = 3.0; Y4 = 9.0; Y5 = 1.5 and a5 = 1.5.

S.B. Welch, S. Salhi / European Journal of Operational Research 102 (1997) 302-319

311

As can be seen in Table 1, the results of Dispfloc are on the whole relatively quick, even for large networks of high density. Here the notation of Min, Mean and Max refer to the minimum, mean and maximum computing times respectively.

4.2. p facility results The accuracy measures used are set at tOldel = 0.05 and tolerr = 0.2 km and the parameters at: /31 = 1.2; /32 = 6; /33 = 1.2; /34 = 1.5;/35 = 3.0 and/36 = 2.5, and Yl, Y5 and c~5 are defined as before. To test the effectiveness of Eq. (6) the complete enumeration of the endpoints was replaced with the scanning of some subset at each iteration. Preliminary tests showed a half scanning technique to be the best compromise between solution quality and computational effort. Here, only 3.5% of problems fail to show an improvement in computing time and of those that do, the time is improved by an average of 50.1%. Of the networks tested, only 8.8% show a change in wp, and those that do only change by on average 1.7%. In our analysis this modification is used in the heuristics P-subdfloc, SA1, SA2 and SA3.

4.2.1. Solution quality The results are displayed in Table 2. The observed values are the mean values of wp per facility, rather than the sum normally used in the p-maxisum problem. This is used to enable easier and more consistent comparison between results of differing values of p. The figures in bold signify the largest mean value of wp over all the heuristics tested. The starting solution of each Simulated Annealing algorithm is obtained using the heuristic given in brackets. An average increase in wp of 5.8% is found by using P-subdfloc compared with P-dfloc, with 5.8% for p = 3 and 6.6%, 6.8% and 4.1% for p = 4, 5 and 6 respectively. This increase in the mean of wp implies that the introduction of the p-maxisum objective whilst keeping the p facilities within the bound set by the maximin primary objective is beneficial when minimising risk. The results of the Simulated Annealing algorithm SA! seem to be very similar to those of P-subdfloc. In fact, the 160 cases tested gave only an increase of 1.2% in the mean value of wp. This relatively small average increase seems to vindicate the assumptions made in the algorithm P-subdfloc. In fact this average increase has the value it does, only due to a few large increases. These special cases occur when a large subset of the facilities are to be sited on a small subset of arcs that are close together and relatively long. This implies that a large number of combinations of possible feasible locations exist at each iteration that are difficult to optimise via a constructive heuristic such as Dispfloc. To overcome this, the algorithm may be improved by detecting such special cases and then guiding the search towards a better solution. These results do seem encouraging in their resemblance to those given by P-subdfloc, and this is re-iterated by the results of SA2, which show no significant difference to those of SA1. The post-optimiser SA3 improved on both the results for SA1 and P-subdfloc with an average increase in the mean of wp of 1.7%. Again this improvement is only due to a small number of special cases which mainly occur in networks having a high number of arcs. Using empirical testing, this refinement has shown that our constructive heuristic P-subdfloc is a good performer as the additional improvement in wp may be considered not so important given the extra computing effort required.

4.2.2. Computing effort The mean computing time for p = 0.2 is illustrated in Fig. 2. Here, it can be observed for all cases studied that time increases exponentially as the number of facilities to be sited increases. It can also be seen that the mean computing time of P-subdfloc is about 2.5 times that of P-dfloc. This increase seems large, but does compare favourably with the computing times for the Simulated Annealing algorithms. Similar exponential growth was obtained for p = 0.5, although the times observed were slightly larger.

312

S.B. Welch, S. Salhi / European Journal of Operational Research 102 (1997) 302-319

KEY timo(secs)

............... ........ ...... .......

2000-

P-dfloc P-subdfloc SA1 SA2 SA3

/ /

!

Table 1 Computational time for Dispfloc (secs)

I

I

i

1500

i

•

i

i i

i

."

i

500

~ ~I o

i

2

p

Min

Mean

Max

20 to 40

0.2 0.5 0.8

8.44 9.84 14.40

18.76 18.09 22.86

25.71 34.25 42.63

60 to 80

0.2 0.5 0.8

26.30 19.60 22.41

28.78 38.60 47.10

62.30 81.90 91.57

100 to 120

0.2 0.5 0.8

30.19 37.08 44.45

62.29 171.55 131.21

87.40 320.12 187.34

.'

'" 7 i ." / i j. .' • • ./" " / ,,/ ; .'// ....." . .J .. • l' °'s .." .'"

lOOO

No. of nodes

.'

~..-" - - - - ." ............. .... / .<,-,.." ~ .... ' " , , P

I

3

4

5

6

Fig. 2. A graph to illustrate the computing times (p = 0.2).

Table 2 Statistical performance of heuristics using half scanning p=0,2

p=0.5

p

3

4

5

6

3

4

5

6

P-dfloc

Min Mean Max

7.27 10.84 16.16

7.10 10.10 15.08

7.10 9.66 13.01

7.05 9.44 12.99

8.31 11.47 16.23

7.56 10.92 16.11

7.32 10.59 15.30

7.32 10.41 15.23

P-subdfloc

Min Mean Max

7.65 11.43 17.42

7.58 10.87 16.13

7.43 10.34 15.06

7.43 9.78 13.18

8.39 12.14 17.72

7.56 11.40 16.88

7.56 11.25 16.51

7.56 I 1.00 15.61

SA 1 (P-dfloc)

Min Mean Max

7.70 11.43 17.23

7.68 10.95 16.08

7.16 10.44 15.30

7.08 10.05 13.73

8.31 12.10 17.69

7.63 11.54 16.88

7.58 11.32 16.56

7.56 11.12 15.61

SA2 (P-dfloc)

Min Mean Max

7.70 11.43 17.23

7.68 10.95 16.08

7.16 10.51 15.31

7.08 10.06 13.74

8.31 12.12 17.70

7.63 11.54 16.88

7.58 1 1.32 16.56

7.56 11.13 15.61

SA3

Min

7.65

7.68

7.59

7.47

8.39

7.66

7.58

7.58

(P-subdfloc)

Mean

11.48

11.00

10.62

10.06

12.20

11.72

11.43

11.19

Max

17.42

16.13

15.30

13.70

17.77

16.88

16.63

15.64

O v e r a l l , it h a s b e e n o b s e r v e d t h a t o t h e r t h a n in a f e w s p e c i a l c a s e s , t h e c o n s t r u c t i v e h e u r i s t i c P - s u b d f l o c p r o d u c e s g o o d s o l u t i o n s f o r t h e p r o b l e m s t e s t e d w h e n c o m p a r e d to t h o s e f o u n d b y s o m e v a r i a n t s o f S i m u l a t e d A n n e a l i n g , b e s i d e s b e i n g r e l a t i v e l y faster.

S.B. Welch, S. Salhi/European Journal of Operational Research 102 (1997) 302-319

313

5. Conclusions and suggestions The investigation developed a heuristic procedure for the single obnoxious facility case. Initially, this involved the building up of circles centred at the nodes of the network and then the use of polygons enabling the inclusion of a pollution dispersion model. The p facility maximin case was then introduced, where the facilities to be sited interact in some way. The heuristic employed a reduction of the continuous feasible location set at each iteration, to that of a discrete set of points. These points were then given a priority rating. This formulation was extended by setting the optimal value for the p-maximin problem as a constraint to a p-maxisum formulation. The heuristic was then tested against two recent versions of Simulated Annealing algorithms. It was found that the Simulated Annealing results showed no significant difference to the results obtained previously, thus reinforcing the claim that our heuristic yielded good approximations to the optimal solutions of the problems tested. In this study we used a specific and well known dispersion equation, but the proposed method can be easily adapted for other models. The extension of the p facility maximin network location model to the plane is worth investigating. The problem may be constrained by imposing restrictions on some forbidden areas, etc. Economical constraints related to accessibility may also need to be considered. The single facility maximin formulation which uses a fathoming technique similar to that highlighted by Karkazis and Papaditriou [14] may be useful as a basis for investigation. The authors are currently investigating this issue and preliminary work has showed promising results.

Acknowledgements The authors would like to thank the referees, whose comments improved the presentation of this paper.

Appendix A. The dispersion equation of Karkazis and Papadimitriou A polar co-ordinate system will be introduced to overcome the rotation used to align the wind direction along the x-axis.

C ( r , h , qb, O) = ~ru( O)q(r,Q¢b)v(r, qb) *exp

(l((rsin~b) -

2

~ q(r, ~b)2 + v(r, dp)-------

,

0 E [0, 27r),

(A.I)

where

q(r, ~b) =0.0804 * (r c o s ~ ) 0'8929,

(A.2)

c,(r, ~ ) = 8.3913 - 1.069 * ( r c o s ~b) 1/2 + 19.445. (r cos ~b) I/3 - 27.3357 * (r cos q~) 1/4,

(A.3)

where: C = The pollutant concentration. Q = Pollutant emission rate (g/sec). u = Wind speed (knots). h = Equivalent stack height. Defined as the sum of the stack height and the rising of the effluent gases due to heat. Valid only when h ~< 300 metres. The direction of the wind being considered subtended at the x-axis. ,~= The angle subtended at 0 by the direction being considered. F = Distance (metres) from the facility. Valid only when r C [ 100, 120000] (metres). This is known as the active interval. The functions q(r, qb) and v(r, dp) were derived to specified accuracy limits; see Karkazis and Papadimitriou [ 14] for details.

S.B. Welch, S. Salhi/European Journal ~f OperationalResearch 102 (1997) 302-319

314

It should be noted here that when the problem is formulated considering the pollution dispersion equation

C(r, h, qb, 0), the topic of what happens to the function as ~b changes must be addressed. Given some direction 0j, then the radius at which some pollution concentration C occurs can be obtained when 0 = 01 and ~b = 0 or 0 = 01 i ~1 and ~b = ~bl for small values of ~bl. It therefore needs to be examined for what values of ~b do the largest and smallest radii occur (depending on which side of the root rh is being studied). In fact it can be proven that except for very large changes in wind speed, the change due to ~b can be ignored; see Appendix B.

Appendix B. The dispersion equation with respect to ~b The following proof assumes that only those values of r outside rh will be used; if the two polygon technique is used, a similar proof can be obtained for values inside rh.

B.1. Proof" the change of C due to ~b can be ignored Proof. Since u(O) is approximated by the straight line side of a polygon, it follows that it has the form y = ax + b. Using polar co-ordinates this gives

u(O) = b / ( s i n O - a c o s O ) ,

0E [0,2~').

(B.1)

Now Eqs. ( A . 1 ) - ( A . 3 ) give

ac

C(r,h,~,O)

O0-

b

u(O) (cos 0 + a sin 0),

OC 3_~=C(r,h,c~,O)[ (

-1

+

(B.2)

y ( r , ~b)2"~

.

Oq

+

( v_~, 1

+

h2

) . Ov _ y(r, qb)

q(T V

Oy

]

(13.3)

It can be assumed that the original problem has been rotated so that the direction being considered is as given below. Now the aim is to compare (Oc/30) when 0 = ~b and ~ = 0, with (Oc/3dp) when 0 = 0 and ~b = ~b, where ~b ~ 0. Also, without loss of generality it can be assumed that ~b > 0. Substituting these into Eqs. (B.2) and (B.3) gives --

30

-

-

3C c?c~

~

,~

¢r(_-~)q(r,O)v(r,O) Q zr(b)q(r,O)v(r,O) x

*exp

-

. exp ( _ ~

dq(r, p--~+q(~,&)3J

v(r-~-0) 2

* - - , -a

( y ( r , ~ b ) 2 ) ) . exp ( _ ~

\q(r,O) z

"3-~+

qb-----f+v(r,(b) v(r, -----~

(B.4)

he

(v(r~)2))

"O-~- ga)-----'-'-~'~ q(r,

"

(B.5)

Consider two cases: a > 0 and a < 0. There is no need to consider a = 0, as considering the rotation that has been undertaken, this would imply that the polygon side is parallel to the direction being considered, which is unrealistic. aO: If a > 0, then (3C/30) < 0, implying that r will decrease. Therefore, it is sufficient to prove tOC/30[ < [3C/O~bI and since both are negative, this reduces to proving (0C/30) > (OC/Odp). To prove this, a proof by contradiction will be used. Therefore, assume (0C/30) < (0C/3(~) is true. This therefore gives the following, since C(r, h, 0, 0) > 0:

S.B. Welch, S. Salhi /European Journal of Operational Research 102 (1997) 302-319

315

l+a~b [ ( y(r, qg)2"~ ( - 1 h2 ) Ov r2q9 ] - - - 2 - - > ~ 0.89295 -1 + -q-~,~T J v(-2,~) + -v(r, - - - -cb) 5 " - ~ + q(r--,-~) 2 J a~b (-0.8929 + 138.1309q~2r °2142 + 154.6991r °2'42) + a

[(

(B.6)

--1 h 2 ) O~] v~,-qb) + -------T v(r, oh) " ,

(B.7)

where, since ~b ,-~ 0 and r E [rh, 120000], we have 0.90 < a < 1.0. Now let the notation (~) imply that the minimum of the replaced function is a and the maximum b. The right hand side then becomes ~0

- - + - cxqg('i38.1309~2r°Xle2+154.6991r 0"2142) + ffq~[(v(~lq~) v(r,q~)3). (

-0 9 91 . (B.8)

By recognising the fact that h E [0, 300] and r C [rh, 120000], where rh is dependent on the value of h, it can be shown that/3 ~ 0 (see Section B.2). This therefore reduces Eq. (B.7) to

(1 + ad?) /a > 154.6991a(br °2142.

(B.9)

Now this may be true when a < 1, but when a > 1 the above is false. Therefore, a contradiction has been obtained. Unfortunately, this only holds when a > 1. However, this constraint is not too restricting when considering the rotation of the wind direction implemented, as 0 < a < 1 would imply that the change in the wind speed is extremely high. If such a model was used, then the approximation techniques used within the proof imply that further investigation is required when 0 < a < 1. []

B.2. Proof" fl ~ 0 Proof. From Eq. (B.8) and v(r,q~) ~ v(r, 0)~ Y A

" -1 h2 /3~a
"

-19.1126"~_0.8929] -3.1297 ] _"

(B.10)

Now y is at a minimum when h = 0 and r = 100, since v(r,O) is monotonically increasing over [ 100, 120000] and h C [0, 300]. Therefore, using the same notation as before, a = -0.6608. This gives a lower bound over h E [0, 300] and r E [rh, 120000]. An upper bound can be found by observing different values of h and setting v(r, 0) to be as small as possible. As only those radii outside rh are being considered, it follows that this occurs when r = rh (v(r,O) is monotonically increasing). However, at the root rh, Karkazis and Papadimitriou [ 14] state

v'(r,O) -v(r,O) - .

(

h2 v(~, b) 2

) 1

q'(r'0) = 0 q(r,O~

(B.11

giving h2 v(r, 0) 3

1 v(r, 0) -

-

-

0.8929 r. v'(r,O)"

(B.12)

Eq. (B.12) is at a maximum when r.v'(r, 0) is at a minimum. This occurs at r = 100, giving r.v'(r, 0) = 3.1295. Substituting this back into Eq. (B.10) gives /3 ~ 0. []

316

S.B. Welch, S. Salhi /European Journal of Operational Research 102 (1997) 302-319

Appendix C. The algorithms In this A p p e n d i x we present the four algorithms: Dispfloc, Facloc, P-dfloc and P-subdfloc.

Dispfloc Step 1. For each node Ni, input a Cartesian co-ordinate and a weighting (xi, Yi, wi). Input arcs (Ni, Nj) G A such that xi < xj. Define the error bounding, tol, for the value of the weighted pollution level, wp.

Step 2. Input the original polygon Pi, l, l = I . . . . . 8, defining the wind speed. A clamped cubic spline technique is used to form a wind speed rose.

Step 3. An initial value for wp is found using a combination of an initially input value and a scanning technique. Set deltastep = wp/yj

Step 4. Calculate the maximum perpendicular error for wind speed, 6w, and an upper bound on radial wind speed error, ew Step 5. For each node Ni C G, a polygon PiJ whose vertices are defined by a wind speed and a pollution level, PL = w i / ( w p * 1000.0), i = 1. . . . . m, is drawn centred at (xi, Yi). Determine the maximum perpendicular distance error 6d and the maximum radial distance error ca. For each nodeNi, iG {1 . . . . . m} and for each (j, k) GA; If ifeasp(arc)= 1 then skip to next arc. Find intersection points ( x , y ) G Pi.I N (j, k) by solving the relevant equation. Store co-ordinates along with the originator node N i and arc (j, k). Any node Nj G pconvN~(Pij) is stored in a record for node Nj. Step 6. Order the intersection points in terms of increasing x, for each arc. Step 7. Along each arc ( i , j ) G A, use a binary array B of length m to check for feasibility. For each node, Ni, i C {1 . . . . . m}; for each ( i , j ) G A For each polygon PkJ; If Ni C pconvNk(Pk,I) then let Bk = 1 (k-th entry of B) for arc ( i , j ) . End polygon; For each intersection point on ( i , j ) from node k; Set Bk ~ Bk + I (mod 2) If ~-~,,z Bn = 0 then the point is noted, as a feasible interval FA has been found. n= 1 The next loop closes the feasible interval. All infeasible arcs ifeasp(arc) are set to equal 1. Step 8. Scale e o according to the polygon's intersections with the feasible arcs. Step 9. If any feasible intervals FA are found then compute final using Eq. (6). If ej < .final < tol then STOP and the solution set has been found. Otherwise: If final ~ e o then: If final < tol and e d < tol then: If this is the 3rd time this has occurred in a row then set wp = wp - deltastep and GOTO Step 5 else: set deltastep = deltastep/y2; wp = wp - deltastep and GOTO Step 5 For each side j of the polygon Pi.t: if 6 ~ > 6w/a5 then a vertex is added to give Pi.l+l. If this is the first time this has occurred since previous infeasible solution then: set deltastep = deltastep/y2 Set wp = wp and GOTO Step 4 If ed < final ~< 2 , ed then: For each side j of the polygon PiJ: if &vj > ~w/a5 then a vertex is added to give Pi.t+l. If this is first time this has occurred since previous infeasible solution then: set deltastep = deltastep/y3 ; Set wp = wp + deltastep and GOTO Step 4 Otherwise: Set wp = wp + deltastep and GOTO Step 5 If no feasible intervals are found then: set wp = wp - deltastep and deltastep = deltastep/y5 If this is the 3rd consecutive infeasible solution then set deltastep = deltastep * Y4 GOTO Step 5 END

S.B. Welch, S. Salhi/European Journal of Operational Research 102 (1997) 302-319

317

Some comments: Note wp is multiplied by 1000.0 within Dispfloc to change the scaling from kilometres to metres, necessary for the dispersion equation and final is as stated in Eq. (C.I):

maxi,j,k final =

(i,j)

(i,j)

d(FAk,1 ,FAk,2 )2 2.0

)

,

( i , j ) C A,

k = 1. . . . . IFA~'j)I.

(C.1)

P-dfloc Step Step Step Step Step Step Step

1. 2. 3. 4. 5. 6. 7.

As Dispfloc. Define tol for deltastep, tolde j and for the error, tolen-. Input the weightings vj, j = 1. . . . . m As Step 2 of Dispfloc. As Step 3 of Dispfloc. Set deltastep = wp/y, Calculate the bounds for the wind speed perpendicular errors 6wa, 8wb and radial errors Ew,,, ewb. As Step 5 of Dispfloc except u(O) + Eb~ is used for the wind speed rather than u(O) As Step 6 of Dispfloc. As Step 7 of Dispfloc. If no feasible intervals FA are found on any arcs then: Set wp = wp - deltastep and deltastep = deltastep/fll If this is the 3rd time in a row that an infeasible solution has occurred then set deltastep = deltastep * f12 GOTO Step 5 Step 8. Each feasible solution interval is reduced with respect to Ed,, e0h, the relevant nodes weighting and the angle subtended by the intersecting polygon and the relevant arc. The maximum error em is then found. If no feasible intervals FA are found then: Set wp = wp - deltastep and deltastep = deltastep/y5 If this is 3rd time in a row that an infeasible solution has occurred then set deltastep = deltastep • [32 GOTO Step 5 Step 9. See Facloc.

Step 10. If no feasible p facility location is found then: For each side j of the polygon Pi,t: if max{SJw,, 6JWb} > max{Sw,, 8Wb}/aS then a vertex is added to give Pi.t+l. Set deltastep = deltastep/fl3, wp = wp - deltastep and GOTO Step 4 Otherwise: If deltastep < toldcl, e m < tole~ and the previous iteration was infeasible then the solution set has been found. If deltastep < tOldeb e m < tole~v and the previous iteration was feasible then wp = wp + deltastep, GOTO Step 5 If deltastep > 4 * tOldel then: If em > tole~v then: If previous iteration was infeasible then deltastep = deltastep/~4 More vertices are introduced to each polygon PiJ to give Pi,t+l. Set wp = wp + deltastep and GOTO Step 4 else set deltastep = deltastep/[35, wp = wp + deltastep and GOTO Step 5 If toldel < deltastep ~< 4 * tOldeI then: If em > tole~T then: If previous iteration was infeasible then deltastep = deltastep/fl6 More vertices are introduced to each polygon Pi,I to give Pi,l+l Set wp = wp + deltastep and GOTO Step 4 elseif this is first time this has occurred since previous infeasible solution then set deltastep = deltastep/fl5 Set wp = wp + deltastep and GOTO STEP 5 If deltastep < tOldeI and em > tole~v then: More vertices are introduced to each polygon Pij to give Pi/+l GOTO Step 4 END

318

S.B. Welch, S. Salhi/European Journal of Operational Research 102 (1997) 302-319

Facloc Step 1. Order the facilities in order of increasing weight. Set PI = 1 Step 2. For each endpoint of the feasible interval lengths FApI, a value is assigned subject to Eq. (6). Step 3. Order the endpoints in order of decreasing objective function. Step 4. Site PI facility at next remaining possible feasible endpoint, subject to Eq. (6). Step 5. Find the remaining feasible interval lengths FApI+l subject to an inverse pollution spread about the facilities already sited using vi, i = 1. . . . . PI - 1, and a normal pollution spread using uj, where j = PI. Step 6. (a) If no feasible interval lengths remain then: Unsite PI facility. (b) If no further endpoints are available then: Set P I = P I - 1 If PI = 0 then no feasible solution is found. GOTO Step 6b else GOTO Step 4. else Set P I = P I + 1 If PI = p then a feasible solution is obtained. GOTO Step 2.

END

P-subdfloc As P-dfloc and including: Step la. Set k = 0 Step 9a. Before Step 9, reduce the set of feasible intervals FA subject to those k facilities already optimised. Step lOa. (Replaces equivalent If statement) If deltastep < tol~lei,e m < tolerr and the previous iteration was infeasible then: Site (k + l)st facility subject to: (a) Order facilities in order of increasing weighting (b) The first facility is sited on an arc that no longer occurs in the set FA when wp increases. If no such arc exists then site first of the ordered facilities at its present location. Set deltastep = deltastep * (p - k) a * ~, k = k + 1 Goto Step 5

References [1 ] R. Batta, S.S. Chiu, Optimal obnoxious paths on a network: transportation of hazardous materials, Operations Research 36 (1) (1988) 84-92. [2] R.L. Church, R.S. Garfinkel, Locating an obnoxious facility on a network, Transportation Science 12 (2) (1978) 107-118. [3] B. Dasarathy, L.J. White, A maximin location problem, Operations Research 28 (6) (1980) 1385-1401. 14] K.A. Dowsland, Some experiments with Simulated Annealing techniques for packing problems, European Journal of Operational Research 68 (1993) 389-399. [5] Z. Drezner, G.O. Wesolowsky, A maximin location problem with maximum distance constraints, AIIE Transactions 12 (3) (1980) 249-252. [6] Z. Drezner, G.O. Wesolowsky, Location of multiple obnoxious facilities, Transportation Science 19 (3) (1985) 193-202. [7] M.C. Fonseca, M.E. Captivo, Location of semiobnoxious facilities, Euro XII, Helsinki, June-July, 1992. [ 8] J. Halpem, Finding minimal center-median convex combination (cent-dian) of a graph, Management Science 24 (1978) 535-544. 19] G.Y. Handler, Medi-centers of a tree, The Leon Recanati Graduate School of Business Administration, Tel-Aviv University, Israel, 1976. [ 10] J. Karkazis, T.B. Boffey, Modelling pollution spread, Studies in Locational Analysis 7 (1994) 91-104. [ l 1 ] J. Karkazis, T.B. Boffey, N. Malevris, Location of facilities producing airborne pollution, Journal of the Operational Research Society 43 (4) (1992) 313-320.

S.B. Welch, S. Salhi / European Journal of Operational Research 102 (1997) 302-319

319

I 12 ] J. Karkazis, E Karagiorgis, A method to locate the maximum circle(s) inscribed in a polygon, Belgian Journal of Operations Research 26 (3) (1986) 3-36. 113] J. Karkazis, P. Karagiorgis, The general problem of locating obnoxious facilities in the plane, in: Proceedings of the 1 lth IFORS Conference, Buenos Aires, 1987, pp. 703-717. 114] J. Karkazis, C. Papadimitriou, A branch and bound algorithm for the location of facilities causing atmospheric pollution, European Journal of Operational Research 58 (1992) 363-373. I15] S. Kirkpatrick, C.D. Gellat, M.P. Vechhi, Optimisation by Simulated Annealing, Science 220 (1983) 671-680. [ 16] E. Melachrinoudis, T.P. Cullinane, Locating an undesirable facility within a geographical region using the maximin criterion, Journal of Regional Science 25 (1) (1985) 115-127. [ 171 E. Melachrinoudis, T.P. Cullinane, Locating an undesirable facility with a minimax criterion, European Journal of Operational Research 24 (1986) 239-246. [18] I.D. Moon, S.S. Chaudhry, An analysis of network location problems with distance constraints, Management Science 30 (1984) 290-307. 1191 D. Moon, L. Papayanopoulos, Minimax location of two facilities with minimum separation: interactive graphical solutions, Journal of the Operational Research Society 42 (8) (1991) 685-694. 120] I.H. Osman, N. Christofides, Capacitated clustering problems by hybrid Simulated Annealing and Tabu Search, International Transactions in Operational Research 1 (3) (1994) 317-336. 121 ] E Pasquill, The estimation of the dispersion of windborne material, Meteorological Magazine 90 ( 1961 ) 33-49. [221 C.R. Reeves, Modern Heuristic Techniques for Combinatorial Problems, Wiley, New York, 1993. 123] S.B. Welch, The obnoxious facility location problem, School of Mathematics and Statistics, University of Birmingham, 1995. 1241 S.B. Welch, S. Salhi, The p-obnoxious facility network location problem with facility interaction, Pre-print No. 96/1, School of Mathematics and Statistics, University of Birmingham, 1996. 1251 S.S. Ting, A linear-time algorithm for maxisum facility location on tree networks, Transportation Science 18 (1984) 76-84. 126] D.B. Turner, Workbook of Atmospheric Dispersion Estimates, Publ. No. AP-26, US EPA, Office of Air Programs, 1970.