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 x_ij=1 if demand location i is allocated to supply location j, otherwise it is 0, and x_jj=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, s_ij, of servicing customer demand at locations i=1,2..n from the supply sites, j. Here, s_ij=w_ic_ij where w_i is a non-zero weight value indicating the demand at customer site i, and c_ij is the cost of supplying a unit of demand from supply location j to customer i. Typically c_ij 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 s_ij is not equal to s_ji. If Equation 1: is replaced with min{max{s_ijx_ij}} 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 x_jj=1 and then compute the sum of the cost of servicing customer demand, C_j, 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 x_ij=1 if x_jj=1 and s_ij-λ_i<0, otherwise x_ij=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 n^th iteration, t⁽ⁿ⁾, is determined using an expression of the form:

Selection of suitable values for a⁽¹⁾ and the λ_i⁽¹⁾ is important to ensure convergence. Christofides and Beasley (1982) and Beasley (1985) suggest using λ_i⁽¹⁾=min(s_ij) across all j, for all i not equal to j, and a⁽¹⁾=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.