Simple heuristics

Exact solution is possible of the network version of the p-median problem. In this model there are n customer demand sites, each with demand weights wi and p supply sites to be chosen from m<=n candidates. These p “median” sites will lie at vertices of the network (in fact optimal solutions may exist anywhere on the network but a theorem due to Hakimi shows that at least one optimal solution will lie entirely at network nodes, which simplifies the overall problem in many cases). In principle, for small n, m and p, solutions can be found by enumerating the possible configurations, but this becomes infeasible very quickly as their values increase in size. Heuristic solutions may provide good or even optimal solutions in a reasonable time, but evaluating the quality of some of these approaches is not straightforward. Other heuristics, by their very nature, do provide bounds on the quality of their approach to optimality.

A basic heuristic for the network-based p-median problem, due to Teitz and Bart (1968) and also described in Section 7.2.2, Heuristic and meta-heuristic algorithms, is known as an “interchange” method and proceeds in the following manner. Let V be the set of m candidate vertices, then: (i) randomly select p vertices from V and call this set Q; (ii) for each vertex i in Q and each j not in Q (i.e. in the set V but not in Q) swap i and j and see if the value of the objective function is improved; if so keep this new solution as the new set Q; (iii) iterate step (ii) until no further improvements are found. Note that this algorithm has the potential to produce improved solutions by repeating step (i) with a different random set and/or by examining good alternative unoccupied locations (e.g. within sub-regions of the network) as preferred choices for the interchange process ― see Densham and Rushton (1992) for more details.

Lagrangian relaxation

Because the p-median problem is of such importance and generality we will describe it in greater detail, together with a particularly effective solution procedure for larger problems known by the rather cumbersome name of “Lagrangian relaxation with subgradient optimization” or LR procedures. Developments of this type of approach are described by Christofides and Beasley (1982) and Beasley (1985). Essentially, this procedure is a form of “branch and bound” algorithm, whereby upper and lower bounds are determined for the problem, optionally used in conjunction with a form of tree search (branching) to produce steadily improved results.

As with many optimization problems a standard formulation of the p-median problem can be made as follows:

subject to:

In this formulation xij=1 if demand location i is allocated to supply location j, otherwise it is 0, and xjj=1 if a facility is opened at location j otherwise it is 0. These requirements can be expressed as an additional constraint equation, requiring that

Equation 1: is the objective function, the total (or average) cost, sij, of servicing customer demand at locations i=1,2..n from the supply sites, j. Here, sij=wicij where wi is a non-zero weight value indicating the demand at customer site i, and cij is the cost of supplying a unit of demand from supply location j to customer i. Typically cij is the shortest path distance across the network between i and j. This information may be computed in real time or supplied as a cost/distance matrix or list. Note that in general sij is not equal to sji. If Equation 1: is replaced with min{max{sijxij}} then the formulation is known as the p-center problem (which is also known to be NP-hard). In this latter case we are trying to minimize the maximum cost of servicing demand, for example to meet pre-defined service criteria.

Equation 2: states that every customer is served by exactly one site. An alternative formulation, 2', is shown as it is used later in this subsection.

Equation 3: states that if a facility is not located at site j then customer demand cannot be allocated to it.

Finally Equation 4: states that a total of p facilities need to be allocated.

The p-median problem on a network can be solved exactly for reasonably large problems (e.g. 1000 customer sites and 50 depots) using the LR algorithm. The solution procedure involves solving a simpler problem (a relaxation of the original problem) that does not satisfy all of the constraints (e.g. fails to satisfy Equation 2:) and then using this solution to obtain upper and lower bounds on the optimum solution. These bounds are then systematically narrowed (the subgradient optimization phase) until, ideally, the solution upper bound (which does satisfy the constraints) meets the lower bound (exactly or to within a given tolerance), providing the optimal solution.

The simpler (relaxed) problem to solve actually looks more complicated. It involves adding a set of n variables (Lagrangian multipliers, λi>0) into the original Equation 1:, “relaxing” the constraint in Equation 2:, as follows:

Observe that the last term in brackets in this expression is the constraint D, which now forms part of the objective function. The problem is ‘relaxed’ because the constraint in Equation 2: has been removed, so it is possible for a customer to be served by/assigned to more than one site — and indeed, when one comes to solve the relaxed problem this frequently occurs. Also note that the final summation can be expanded and combined with the initial summation, as follows:

Minimization of the main term in brackets above is thus a simple enumeration task using the first term of the expression, since the second summation is a constant. We select each of the potential sites, j, in turn by setting xjj=1 and then compute the sum of the cost of servicing customer demand, Cj, using the selected facility minus the Lagrangian, λi, i.e.

If this sum is greater than 0 it is set to 0 otherwise the value is stored and the computation proceeds for the next j. The p smallest (largest negative) values obtained in this way are taken as the locations of the initial set of p-median sites. The allocation of customer sites to these selected locations is made by setting xij=1 if xjj=1 and sij-λi<0, otherwise xij=0. This allocation will not generally satisfy the original constraints, but the computed total cost obtained using Equation 1: will provide a lower bound, LB, to the true optimum solution.

We now obtain an upper bound, UB, by allocating the customer sites to their nearest supply location and again computing the total cost using Equation 1:.

To obtain an optimal or near-optimal solution the values of the initial Lagrangian multipliers need to be adjusted (reduced) in such a way as to eventually satisfy the original constraints. The updated values for these multipliers are obtained using a step adjustment to their values during the previous iteration:

The stepsize value at the nth iteration, t(n), is determined using an expression of the form:

Selection of suitable values for a(1) and the λi(1) is important to ensure convergence. Christofides and Beasley (1982) and Beasley (1985) suggest using λi(1)=min(sij) across all j, for all i not equal to j, and a(1)=2, with procedures for reducing this in subsequent iterations (e.g. progressive halving) depending upon the convergence behavior. Note that the term in the denominator of this expression is the sum of squared ‘violations’ of the constraint we have relaxed, so in this case it is based on the number of customers who have been assigned to more than one supply center.

After each iteration new UB and LB value are obtained as before and the process continues until satisfactory or complete convergence has been achieved. An optional standard tree search procedure can be included at this stage if the process fails to converge. Daskin (1995, pp 225-228) provides a fuller explanation and simple numerical worked example of this procedure. For large problems the overhead of computing the initial upper bound has been found to be high, and as an alternative a simple heuristic such as the interchange method described earlier can be used. Indeed, as we have noted, the heuristic itself can sometimes produce optimal solutions directly.

Performance of the LR algorithm for a range of problems is described by Beasley (1985) and for heuristic implementations of essentially the same procedure, by Beasley (1993). A set of test data for the problems described by Beasley in these articles and a range of related problems (e.g. other problems involving depot and warehouse location, with or without capacity constraints) is available at his web site (see Appendix for details). TSPLIB can also be used as a test dataset for p-median software.

Problems that involve <500 vertices (customers) and modest values of p (e.g. p<50) can be solved directly and optimally by commercial mixed integer program (MIP) solvers such as CPLEX and LINDO. Larger problems remain difficult for such software and/or may require an excessive amount of time to complete — typically an hour or more of run time. Beltran et al. (2004) use the TSPLIB dataset in their semi-Lagrangian solution approach which is built on similar ideas to those of Beasley. Their tests were implemented in MATLab in combination with the CPLEX library V8.1. They use this approach to solve large p-median problems (3000+ customers and up to p=500 supply locations) to optimality or within 1% of optimality. More recently Avella et al. (2006) have reported results on their earlier study of their “Branch and Cut and Price” (BCP) algorithm. BCP can be used to solve large p-median and related problems in conjunction with CPLEX or similar library software. They have tested their methodology against the TSPLIB and OR-Library datasets, with very encouraging results, greatly outperforming CPLEX used on its own.

Comparison of alternative p-median heuristics

The S-Distance software project, mentioned earlier in this subsection and described further in Sirigos and Photis (2005), provides tools to compute solutions to quite large p-median and p-center problems using a range of heuristics, together with solutions for selected coverage problems. The heuristics they implement include many of those described in Taillard (2003), Daskin (1995) and Dyer and Frieze (1985). Data may be read from various file formats, including OR-Library files and network node/link files in dbf format. These various methods may then be compared in order to obtain an estimate of their relative efficiency. For the purposes of this exercise we used an S-Distance test dataset comprising the street network in Tripolis, in the region of Arcadia, Southern Greece. The data consist of 1358 vertices and 2256 edges. Each vertex or node has an integer demand value associated with it, ranging from 1 to 137, plus 35 nodes with 0 demand. Each edge or street link has a “to-from” and/or “from-to” cost of travel assigned to it, with selected edges being designated as one-way. In this test we sought 5 locations to service the demand using several different heuristics.

Figure 7‑15A to Figure 7‑15C show the results of running these procedures on the test dataset. In these figures the red circles identify the solution points (median locations) and the darker lines show the network edges that each center uses to service the demand. The figure in brackets in each case shows the objective function value, z, achieved, in millions of units (sum of demand times cost of travel). The basic greedy-random heuristic used by S-Distance is that defined by Dyer and Frieze (1985). This simply involves selecting the highest weighted demand point (or any at random if there are multiple locations with equal weights). The next location is then chosen by selecting the demand point that has the largest weighted distance from its nearest center. The algorithm is of complexity O(np) and yields results that are within 2-3 times the optimum (and so are not very good!). By contrast, the CLS procedure often produces very good results (within 1% generally) but requires longer to run, with a complexity of at least O(np(n+1)). CLS is essentially a form of alternating location-allocation (ALT) procedure, like that of Cooper described earlier. Initial solutions are then perturbed according to a greedy interchange rule that attempts to swap each ALT solution point in turn with another candidate site, and examines whether the resultant solution improves the ALT local optimum.

The greedy heuristics shown in Figure 7‑15A and B ran in almost no processor time, with the greedy-add algorithm producing a much improved z-value. The five locations it chose, however, are very different from those in Figure 7‑15C, which was the solution found by three other algorithms: candidate list search (CLS), Lagrangian relaxation (upper bound solution, LR), and variable neighborhood search (VNS). Of these VNS was fastest to run, requiring under a second of processor time, whilst CLS took around 3 seconds and LR over 3 minutes. However, the LR procedure provides both upper and lower bounds on the solution, which in this case showed that the lower bound was 1.177, but this solution is infeasible, so it is extremely likely that the solution shown is, in fact, the true optimum. The final image, Figure 7‑15D, shows the allocation of customer demand to the individual centers. The size of the colored circles reflects the demand weights assigned to the network vertices.

Figure 7‑15 Comparison of heuristic p-median solutions, Tripolis, Greece

A. Greedy random (z=1.526) |
B. Greedy add (z=1.190) |

C. LR, VNS and Candidate list search (z=1.180) |
D. Candidate list search, Demand allocation |

Based on data provided by S Sirigos, Dept of Planning & Regional Development, U. of Thessaly, Greece

Similar heuristics may be applied to the p-center problem. The total cost of servicing demand in this case will always be greater than or equal to that for the p-median solution, but the maximum distance (cost of travel) for customers travelling to or from service centers will be minimized. The diagram illustrated at the start of this Guide is the p-center solution equivalent to that shown in Figure 7‑15D. As can be seen, the selected centers are more evenly spread across the City as the algorithm seeks to avoid any very long links between customers and service centers. Particularly noticeable are the locations of the most easterly and south-westerly service centers in this solution. It is interesting to note in this case that this p-center solution is very similar to the p-median solution illustrated in Figure 7‑15B.

These network-based p-median and p-center models can be compared with the discrete location planar equivalent, where demand is located at nodes as before, but facilities may be located anywhere in the plane and the objective function is for median location. In the example shown in Figure 7‑16 the facilities located are approximately optimal with respect to demand using the Euclidean metric rather than network distance. The locations are similar to those identified by the p-center optimization process and significantly different from the best network-based p-median solutions.

Figure 7‑16 Facility location in Tripolis, Greece, planar model

A. Customer demand (nodes) |
B. Median locations and assignments |