Speeding up the optimal method of Drezner for the p -centre problem in the plane (cid:2)

This paper revisits an early but interesting optimal algorithm ﬁrst proposed by Drezner to solve the continuous p -centre problem. The original algorithm is reexamined and eﬃcient neighbourhood reductions which are mathematically supported are proposed to improve its overall computational performance. The revised algorithm yields a considerably high reduction in computational time reaching, in some cases, a decrease of 96%. This new algorithm is now able to ﬁnd proven optimal solutions for large data sets with over 1300 demand points and various values of p for the ﬁrst time. © 2016 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ).


Introduction
The p -centre problem seeks to minimise the maximum distance or travel time whilst ensuring all the n demand points are covered by at least one of the p chosen facilities.This problem can be categorised as either the vertex p -centre problem or the absolute p -centre problem.In the former, which is the discrete case, the optimal facilities are part of a set of the potential facility sites which can be either the demand points or other known sites.However, in the latter the facilities can be located anywhere along network edges (as introduced but not solved by Hakimi (1965) ) or in the plane.
In this paper, we will explore the absolute p -centre problem in the plane, which is also known as the continuous or the planar p -centre problem.It is worth noting that the continuous p -centre problem, besides being used for interesting real life location applications that will be briefly mentioned next, could also provide a greenfield solution which can be used as a guide to identify potential sites for the discrete case as in some cases this data can be very expensive to gather.In addition, the p -centre problem can also be used as a basis for academic research in the general area of global optimisation including other continuous related location problems.
Here are some of the papers describing real-life problems tackled by p -centre models.One of the earliest applications is by Richard, Beguin, and Peeters 's (1990) who used the p -centre problem to locate fifteen fire stations in the Belgian province of Luxembourg.Pacheco and Casado (2004) located a number of health resources such as geriatric and diabetic health care clinics in the rural area of Burgos in Spain.Wei, Murray., and Xiao (2006) adapted their Voronoi-based algorithm developed for the constrained continuous p -centre problem to locate twenty-five emergency warning sirens in Dublin, Ohio.Kavah and Nasr (2011) modified a harmony search heuristic to locate bicycle stations in Isfahan, Iran by solving the conditional and unconditional discrete p -centre problem.Finally, Lu (2013) used the p -centre problem to locate a number of urgent relief distribution centres after the 7.3 Richter scale earthquake in Taiwan.
Most of the real-life applications for the p -centre problem have been solved successfully using powerful heuristics and metaheuristics.However, recent developments in exact methods, with the advances in computing power, memory management and powerful commercial optimisation software such as IBM ILOG CPLEX, mean that the proven optimal solution can now be worth exploring for larger problems.This study aims to respond to such scientific and technological change.In addition, if an optimal solution can be found in a reasonable amount of time, this will provide flexibility in performing scenario analysis for strategic planning purposes which is of extreme importance in practice due to the massive investment usually required.

A brief literature review
The single facility minimax location problem (1-centre) in the continuous space has a long history and was posed originally in 1857 by the English mathematician James Joseph Sylvester .A few years later, in 1860, he proposed an algorithm to solve it.Elzinga and Hearn (1972) developed an efficient and widely used geometrical-based algorithm to solve the problem optimally.Their algorithm was adapted and enhanced by many authors including Elshaikh, Salhi, and Nagy (2015) .For more information on continuous centre problems and references therein, see Drezner (2011) .The idea was extended to find solutions to multi-facility location problems including the p -centre problem.Hakimi (1965) was the first to formulate the continuous 1-centre problem in a network, and Minieka (1970) studied the case where p > 1.The first paper discussing the p -centre problem in the plane was by Chen (1983) .The problem has been shown to be NP-hard when p is variable, see Megiddo and Supowit (1984) .For a fixed value of p the problem can be solved in polynomial time, O (n 2 p+4 ) , though requiring an excessive amount of computational effort especially for larger values of p , see Drezner (1984) .
There exist a few variations of the continuous p -centre problem.For example, Chen and Handler (1993) proposed an efficient algorithm to solve the conditional p -centre problem.Here, the aim is to locate p facilities given that q facilities already exist.Wei et al. (2006) suggested a Voronoi-based algorithm to solve the constrained continuous p -centre problem where the facilities cannot be located within some forbidden regions such as rivers, lakes, military areas etc. Chen and Chen (2013) used Minieka's algorithm and the relaxation method to solve the α-neighbour p -centre problem.In this variation, each demand point is covered by at least α facilities which can be important in the case of facility disruption.
Among the most recent theoretical work is the use of the relaxation concept, where a large problem is broken down into relatively much smaller and more manageable sub-problems that are easier to solve.For more details on this particular topic, see Chen and Handler (1987) , Chen and Chen (2009) and Chen and Chen (2010) .For the discrete case, though not directly related to our research, the following studies by Elloumi, Labbe, and Pochet (2004) , Brandenberg and Roth (2009) and Caruso, Colorni, and Aloi (2003) can be found to be interesting and also informative.In both the discrete and the continuous problems, Cooper 's (1964) Multi-Start method, which is based on the locate-allocate principle, is often used to produce an upper bound for optimal methods or initial solutions for metaheuristics.This paper will be analysing the original continuous p -centre problem by revisiting an interesting, though originally very slow, optimal algorithm proposed thirty years ago by Drezner (1984) .This older method used a subset of facility locations based on specific circles rather than demand points.As this algorithm is the basis of our research, it is detailed in the next section.
The contributions of this study include: (i) revisiting an early but slow optimal algorithm for the continuous p -centre problem; (ii) introducing neighbourhood reduction schemes supported mathematically to improve drastically the computational performance of this exact method; (iii) embedding an adaptive CPLEX policy to further enhance its efficiency; (iv) solving optimally for the first time relatively much larger problems with up to 1300 demand points and up to 100 facilities.
The paper is organised as follows: the investigated exact method is introduced and described in Section 2 , alongside initial results based on the original algorithm.Section 3 proposes the suggested enhancements to the algorithm which are supported by new lemmas and proofs.The computational results are given in Section 4 followed by an adaptive CPLEX policy in Section 5 making this revised optimal algorithm even more efficient.The overall computational results are given in Section 6 .Our conclusions and suggestions are summarised in the final section.

Introduction
Drezner's algorithm is based on the idea of Z-maximal circles .A circle is defined as maximal based on a given upper bound, Z .The set of maximal circles based on Z is then identified and their respective centres are then used as a subset for the potential facility locations.
Let us define the following notations.
I : set of demand points indexed by i = 1 . . .n ; J : set of all possible circles indexed by j = 1 . . .m ; C j : circle j defined by its centre (x c j , y c j ) and radius r j , j ∈ J ; K : subset of I ; R ( K ): the radius of the smallest circle encompassing all points in K ; d i , j : Euclidean distance from demand point i to the centre of circle C j , i ∈ I , j ∈ J ; p : number of facilities to locate; d i,l : Euclidean distance from demand point i to demand point l ; Z : the upper bound at a given iteration; J Z : set of Z-maximal circles ( J Z ⊂ J ).
Definition 2.1.The closure of circle C j is the set of demand points encompassed by circle C j which is defined as Definition 2.2.The minimum covering circle ( MCC ) of the set K is the smallest circle encompassing all points in K with radius R ( K ).
We can now define a Z -maximal circle in the following way, as given by Drezner (1984) .
Definition 2.3.A circle C j with radius r j is said to be Z-maximal (often simply called maximal) if: Drezner proposed two ways to solve the p -centre problem using Z-maximal circles.The first, which will be referred to as CP (a )  0 , uses the set covering problem to find the minimum number of Zmaximal circles needed.First, let the input A i , j be defined as (2) where The objective function (1) refers to minimising the number of Z-maximal circles.Constraint (2) guarantees that every demand point is encompassed, or covered, by at least one Z -maximal circle.( Drezner, 1984 ).In the second method, referred to as CP (b)  0 , a new constraint (4) is added to CP (a )   0 to impose that the number of covering circles has to be equal to p , while the objective function (1) is omitted turning the problem into a feasibility problem.
( CP (b)   0 ): Find x j ∈ {0, 1} ∀ j ∈ J Z subject to (2) and (3) , j∈ J Z x j = p. (4) If the minimum number of covering circles found in (1) is ≤ p or if CP (b)   0 is feasible, then the upper bound is decreased by setting Z to the radius of the largest Z -maximal circle from the obtained solution, and the process of identifying the Z -maximal circles is then repeated.Otherwise (i.e. the minimum number is > p or CP (b)   0 is infeasible), the current upper bound Z is taken as the optimal solution and the algorithm terminates.
Before we use Drezner's optimal algorithm, as described in Fig. 1 , we shall first define the following additional notations.
C 1 J : the set of null circles created from one critical point only (i.e., note: J : the set of circles created from two critical points defining its diameter; C 3 J : the set of circles made up from three critical points forming an acute triangle. It is important to note that an appropriate heuristic must be used to find an initial upper bound in Step 2. For instance, a simple multi start heuristic can be used.In this study we opted for the H 2 heuristic proposed by Drezner (1984) for consistency reasons.

Initial results & the need for an improved implementation
Our initial results were found for two TSP-Library (2015) data sets, namely pr 439 and rat 575 which represent a 439-city problem and a 575-rattled grid problem, respectively.Note that the basic tricks of using squared distances were also adopted here when required to improve code efficiency (e.g. when distances are compared, or for non-acute triangle detection).
Both CP (a )   0 and CP (b)   0 were initially used to solve the p -centre problem, and both were found to take a considerable amount of computational time as a large number of iterations was required.As an illustrative example, we show the result found for the TSP-Library data sets pr 439 and rat 575 where p = 90 .For the data set pr 439, the 90 -centre problem was optimally solved using CP (a )   0 requiring more than 38 hours (i.e.137692.6 seconds) and 4580 iterations.When using CP (b)  0 , the time was reduced to just below 3 hours (10654.30seconds) while using 393 iterations only.For rat 575, an optimal result was obtained using CP (b)  0 , however it required nearly 30 hours (107916.0seconds) and 2729 iterations.Furthermore, when using CP (a )  0 , the program was stopped after the time limit of 2 days with only one feasible solution found with a value of 21.471 (a percentage difference of 18% from the optimal solution).It will be shown later that the optimal solution can be found in less than half an hour (996.43 seconds) with our improved method.This example highlights the importance of developing ways to enhance the efficiency of Drezner's algorithm optimal algorithm.

Modification of the covering problem (enhancement zero)
Traditionally, the continuous p -centre problem is formulated as the Euclidean unweighted p -centre problem.This multiple facility location problem has been examined by a small number of authors, see Plastria (2002) and the references therein.It can also be formulated as a non-linear mathematical programming formulation.However, the formulation that we will use in this paper is similar to Drezner's CP (b)   0 formulation with two commonly used additions consisting of (a) an objective function that aims to minimise the largest radius ( 5) and (b) an extra constraint to deal with the characteristics of the p -centre (8) .This formulation, referred to as CP 1 , will be used throughout this work.
where D : the maximum distance between a facility and a demand point.
The use of CP 1 is first tested on the previous two data sets.It was observed that for p = 90 and 100, the computational times were 1258 and 462 seconds respectively, approximately 9 (resp.7) times faster than using C P (a )   0 (resp.C P (b)  0 ).It is also worth noting that CP 1 has an advantage over Drezner's original suggestions as the optimal solution value D is much tighter leading to requiring a relatively smaller number of iterations.Although it may be harder to solve CP 1 than C P (a )   0 or C P (b)  0 , the last two require a large number of iterations, each including a lengthy Z -maximal circles identification step.

Observations
Optimal solutions for p = 10 , 20 . . ., 100 were found using CP 1 instead of CP (a )   0 or CP (b)   0 for the TSP-Library data sets pr 439 and rat 575.The results are given in the Appendix under Tables A.1 , A.2 and Fig. A.1 .Based on these results, it can be observed that there are two areas where enhancements could be introduced in an attempt to shorten the overall computational time.These include: (a) the way the Z -maximal circles are identified from one iteration to the next; (b) a choice of a compromise between the quality of a feasible solution and its corresponding computational time when solving CP 1 (i.e.finding an optimal solution or just a good feasible solution).
This paper will now investigate several ways in which the original algorithm using CP 1 can be efficiently implemented.Sections 3 -5 will cover (a) and Section 6 will deal with (b).
Note that the introduction of CP 1 , instead of using CP (a )   0 or CP (b)  0 , could be considered as our first enhancement due to generating tighter bounds.However, for simplicity and conciseness, the results of CP 1 will be used as our starting point from which we will base our improvements.

Enhancement one: EHA-based implementation
The Elzinga-Hearn algorithm ( EHA ) is used to find the MCC of a set of demand points.As this is repeatedly needed in Step 3 of the FMC algorithm, in order to calculate R ( Cl j ∪ { i }), two ways in which the overall time performance can be enhanced are highlighted.

Early termination
The EHA starts with a circle made from any two selected points and continues to find a covering circle of increasing size until all points are covered.It is important to realise that in the FMC algorithm, the exact centre point and the radius of the MCC are not needed: we simply aim to establish whether or not the radius of the MCC will be larger or smaller than the upper bound Z .If the MCC is smaller than the upper bound, then the EHA will continue until the end as normal.However, it can be terminated early if the circle's radius exceeds Z during the algorithm.This is because at each iteration in the EHA , the new circle's radius is either the same or larger.Therefore, if a circle has a radius ≥ Z at any point in the algorithm there is no need to continue as the final circle (the MCC ) will be even larger.

More informative initial points
Instead of starting the EHA from random points or selecting points using selection rules, such as the ones adapted by other authors including Welzl (1991) and Elshaikh et al. (2015) , we take into account the information we have already found.In other words, the two or the three critical points that defined the circle found at a current iteration are the points that we choose as our initial points for the EHA .This makes the selection deterministic and yields faster results.
This double enhancement, referred to as Enh1, is incorporated only into Step 3 of the FMC algorithm and does not affect the total number of iterations of the algorithm.

Enhancement two: efficient recording of the Z -maximal circles
At each new iteration in Drezner's algorithm, the process of finding the Z -maximal circles begins again from the start irrespective of earlier iterations.However, when examining the first set of results it was observed that many of the same circles were being classified as Z -maximal during successive iterations.
As an example, Table 1 shows the number of Z -maximal circles found at each of the first 10 iterations of the original algorithm for the data set pr 439 with p = 100 .In this example, approximately 17% of the new Z -maximal circles need to be identified at each iteration only, as the other ones have already been found in previous iterations.Therefore, a technique to identify whether a circle is Z -maximal or not in subsequent iterations is worthwhile constructing.
Lemma 1.If circle C j is Z t -maximal at iteration t , then it is also Z t+1maximal for iteration t + 1 if and only if its radius r j < Z t+1 .
Proof.We know at each iteration t , the upper bound Z strictly decreases.Therefore, we can say Z t > Z t+1 .For circle C j to be a Z -maximal circle at iteration t , the following two conditions need to be satisfied: The information denoting whether or not circle C j has been found to be Z -maximal or not can be stored in a binary or logical vector CircMax where This result is incorporated into Steps 2 and 3 of the FMC algorithm to avoid performing redundant calculations.We will refer to this enhancement as Enh2.

Enhancement three: fast identification of some non-Z -maximal circles
This enhancement, which we will refer to as Enh3, aims to quickly identify some non-Z -maximal circles without performing unnecessary calculations.As an example, take circle C j with a centre point (x c j , y c j ) and radius r j < Z .We can now create a new circle C + j centered at (x c j , y c j ) and with radius Z .Therefore, it is clear that Lemma 2. If s ∈ I is not covered by C j (i.e.s ∈ Cl j ) but is strictly covered by C + j , then circle C j is not Z-maximal.
Proof.Let s ∈ I with s ∈ Cl j but strictly covered by C + j .Then the smallest circle, C , containing s and the whole circle C j , contains all the points in Cl j and is strictly contained in C + j .Hence, C 's radius is at least R ( Cl j ∪ { s }) and is strictly less than Z .It follows that R ( Cl j ∪ { s }) < Z , and so C j is not Z -maximal.
Thus a minimum distance, or threshold, of value Z is established.In other words, if there is at least one demand point not covered by circle C j which lies within this distance, then the circle cannot be classified as Z -maximal.
In summary, if we can conclude that circle C j is not Z -maximal.Additionally, a maximum threshold of 2 Z can also be added using Lemma 3 .Lemma 3. Take any demand point s ∈ I not covered by C j .In case Proof.Take s ∈ I with d s , j ≥ 2 Z .Consider the circle C with centre s and radius 2 Z .As r j < Z , the centre of C j is not encompassed by C .Therefore, the circle arc of C j lying within C is strictly less than half the circle.
But the critical points of C j span at least half the circle, and so cannot all lie within C .Therefore, Thus if a point that lies at a distance ≥2 Z from (x c j , y c j ) is added to the set of points encompassed by the circle C j , the MCC that covers all these points would have a radius ≥ Z .Thus, if this information is known, any point in this area does not need to be checked again and hence computational time can be saved without affecting the quality of the solution.In summary, if then we can conclude that circle C j is Z -maximal.These two observations lead to the construction of a checking area for circle C j , say Check j .This is represented by the shaded area in Fig. 3 , and is defined as follows: (12) We can therefore conclude that further calculations must be performed only if the two observations above are not true and Check j = ∅ .We incorporate Enh3 into Steps 2 and 3 of the FMC algorithm.

Enhancement four: identifying non-Z -maximal circles
If circle C j is not maximal, then there must be a demand point i ∈ Cl j such that R ( Cl j ∪ { i }) < Z .If this point is recorded, in the next iteration this demand point can be the first to be checked and hence repeated computations can be discarded.If the MCC of the next iteration is still < Z , then we can deduce that this circle is still not Z -maximal thus saving computational time.If the MCC is ≥ Z , we either continue with calculations and conclude it is now classified as Z -maximal, or we record the next demand point to cause C j to be non-Z -maximal if it exists.In other words, either way will provide us with useful information that can be used in subsequent iterations.
As an example, say at iteration t it takes q j points to find a demand point that determines circle C j as not Z -maximal.This means the next iteration ought to start with the q th j point instead of starting from scratch at the beginning.This saves the computational time it takes to check the previous (q j − 1) points, say Sa v t j .As this scheme is applied to C j where j = 1 , . . ., m , the saving at iteration t could be significant and of the order of m j=1 Sa v t j .Let Start be an integer vector of dimension m .The entry Start j denotes which demand point i should be checked first in the next iteration to see if circle C j is Z -maximal or not.
This enhancement, referred to as Enh4, is incorporated into Steps 2 and 3 of the FMC algorithm.

Individual performances
The enhancements were first analysed separately so that each one's improvement in computational time could be assessed and its impact measured.For illustrative purposes, the computational times for the individual enhancements for the data set pr 439 where p = 70 , 80 , 90 and 100 are first shown in Fig. 4 .This is then followed by combining all the refinements together using a certain order that will be based on the individual enhancement performances.
Fig. 4 suggests that the best enhancement, giving an average decrease in computational time of 84.42%, is Enh3.By providing minimum and maximum thresholds by which the demand points are checked reduces many calculations as many points sit  outside the checking area.Enh4 yields the second best result with an average decrease of 83.26% in computational time.By starting at the last known non-Z -maximal circle all previous demand points can be disregarded, thus avoiding the unnecessary calculations that they incur.Enh1 is the third best at improving the overall computational time, with an average decrease of 50.65%.This enhancement reduces the number of calculations by terminating the EHA algorithm earlier whenever possible.Also, by choosing the current critical points as the initial points, the EHA will have less iterations to find the MCC .Finally, Enh2 improves the computational time the least.This is due to not dealing with the Z -maximal circle calculations directly; it simply minimises how many circles are needed for these calculations.The average improvement of computational time for Enh2 is 26.26%, which is still significant.

Combined performance
The four enhancements are embedded into Drezner's original algorithm that uses formulation CP 1 .These are added in the order of individual performances observed earlier which is as follows: Enh3 -Enh4 -Enh1 -Enh2.To assess the incremental gain of these enhancements we also conduct the following experiment: in the first run we use Enh3, in the second we use Enh3 and Enh4, and in the third Enh3, Enh4 and Enh1 are used.The fourth run consists of the overall algorithm with all the enhancements incorporated as noted earlier.The results are shown in Fig. 5 .
It is clear that the enhancements greatly improve the computational time.The first enhancement reduces the total computational time by an average of 84.49% as noted earlier, and by adding Enh4 this is decreased further to 90.26%.After the addition of Enh1, the average decrease becomes 96.46% and finally with all enhancements added this reaches a massive saving of 96.71%.In other words, just above 3% of computational time is really needed on average, leading to an exciting and strong result.
It is also worth noting that the incremental decrease in computational time is not directly additive as there is a high level of association between their individual contributions.For instance, after gaining 84% with Enh3, one might expect Enh4 to yield 83% of the remaining 16%.This would therefore give a new decrease of approximately 97%.However, it only decreases it to just over 90% (i.e., an extra 5.8% only).

The complete revised optimal algorithm
The revised FMC algorithm is given in Fig. 6 .It is similar to the original FMC algorithm except Steps 2 and 3 in Fig. 2 have been modified accordingly to accommodate the enhancements described in this study.The revised Drezner algorithm is similar to the Drezner's original algorithm stated previously in Fig. 1 , except that in Step 5 the formulation CP 1 is used instead of CP (a )   0 or CP (b)   0 and an extra step (Step 3 shown in Fig. 7 ) has been added to accommodate the enhancements.For completeness, we reproduce the full revised optimal algorithm in Fig. 7 .

Computational results
The proposed algorithm was coded in C + + on a HP Elitebook 8570w with 12GB of memory.The IBM ILOG CPLEX 12.6 console was incorporated into the program using default parameters.
Tables 2 and 3 show the results found for the data sets pr 439 and rat 575.The first column titled p shows the required number of facilities.The initial upper bound value, denoted by Z in column 2, was found from a 10 0 0 iteration runs of the H 2 heuristic described in Drezner (1984) .The next column, titled Z * , shows the optimal solution value, followed by the computational time (in seconds) required for the revised Drezner optimal algorithm to find Z * in the Loop CPU Time column.Note that this result excludes the computational time consumed by the H 2 heuristic.
Other information, such as how many loops (iterations) are needed to get the optimal solution value, the total time spent on computing the Z -maximal circles and the total time spent on computing the result in CPLEX are reported alongside their corresponding percentages in the remaining columns.(Note that these two individual percentages when added are below 100% due to other calculations.) For completeness, we also produced a summary result in Table 4 to show for both instances and for each value of p the new and the old duration including the percentage decrease.It is clear to see that the enhanced method has greatly reduced the computational time for both data sets.As an example, it took just over 4 hours average computational time for the data set pr 439 previously, whereas now the average time is just over 12 minutes leading to a massive average reduction of 96%.Note that these computational times do not include the computational time for the H 2 heuristic.
For the data set rat 575, the computational time has also been reduced.For the smaller values of p (10, 20 and 30), the majority of the time was taken computing the Z -maximal circles leading to a reduction of over 90%.However, for the other values of p the majority of the computational time is taken up solving the problem in CPLEX leading to an overall relatively small though still significant reduction of nearly 50%.This observation led us to face a challenge that will be explored in the next section.
Furthermore, our findings could be compared to the relaxationbased algorithms of Chen and Chen (2009) for the only reported results for the TSP-Library data set pr 439.In this particular instance, our total computational time (inclusive of the computational time required for the H 2 heuristic) is found to be greater than theirs.However, it is also important to note that our optimal algorithm is deterministic and hence relatively more robust, as it      is not sensitive to several factors including the initial subset of demand points or the number of demand points added to the subset at each iteration.

A self-adaptive CPLEX policy
In this section, we investigate how to balance the time spent between computing the Z -maximal circles and the level of the solution quality which we consider to be acceptable when solving CP 1 .However, to guarantee optimality, we need to show at one stage that CP 1 has no feasible solution and hence the final iteration needs to run to the very end.In other words, it is not possible to reduce the computational time by terminating the search earlier in the last run.
Table 5 shows the total time taken in CPLEX compared to the time consumed in the last iteration in CPLEX.Though a relatively considerable amount of time is used in the last iteration accounting for approximately 10-20% of the total computational time, the computational time taken in the previous iterations is nonetheless worth exploring for possible improvement.A compromise feasible solution to save computational time in CPLEX while limiting the total number of iterations of the entire algorithm will be our focus in this section.
There are several ways in which the search can be terminated early in previous runs whilst producing a feasible solution for CP 1 .An example would be to impose a time limit, however this does not always guarantee that a feasible solution will be found within that time and so other options are investigated.
Our study adopts a strategy by which we manipulate the duality gap so that CPLEX terminates earlier with a good feasible, but not necessarily optimal, solution whenever it manages to find at least one.However, the value of the duality gap can be both sensitive and critical which can make our algorithm less robust.The algorithm cannot terminate too early as it could simply increase the number of iterations greatly, and therefore increase the time spent computing the Z -maximal circles.It is therefore important to find a reasonable compromise that we wish to devise.In this study, we propose the following self -learning CPLEX policy which takes into consideration information from previous iterations.
It is worth noting that the following duality gap policy is only implemented when CPLEX finds at least one feasible solution in any run of CPLEX.However, if no feasible solution has been identified in a given run, CPLEX continues until the maximum time limit is reached where the search terminates.Hence the obtained Z value of the previous run is used as the final solution, which obviously cannot be guaranteed to be optimal.

An adaptive CPLEX policy
At iteration t , the moving average for the computational time for calculating Z -maximal circles ( T Max ) and solving the problem in CPLEX ( T CPLEX ) based on the last α iterations is respectively defined as follows.
where A = { T Max , T CPLEX } , and A t is the corresponding time at iteration t .
We define α as In other words, the classical average is used if t < K , otherwise the moving average over half of the past iterations is adopted.In this study, we used K = 6 based on preliminary results.
We use the following scheme based on the performance ratio then the time for computing the Z -maximal circles is much larger than the time spent solving the problem in CPLEX.Therefore, the number of iterations need to be reduced as much as possible, and so we set the duality gap to 0%.(b) However, if then the majority of the computational time is spent solving the problem in CPLEX, and therefore we wish to exit CPLEX sooner with a feasible solution rather than seeking an optimal one, hence we set the duality gap to be 1%.
(c) If ξ has any other value, then the computational times are considered to be more or less similar.In this case, we wish to reach a balance between finding the near optimal solution and leaving CPLEX early, hence we set the duality gap to be 0.5%.In summary, the following conditions related to the duality gap are given.

Duality Gap
This policy, which uses adaptive learning, is less sensitive to the effect of the data's distribution on the computational time and therefore it is very reliable.
The final results for rat 575, that include the results where the CPLEX adaptive policy is incorporated, are found in Table 6 , displayed alongside the total computational time required to optimally solve this data set using the enhanced algorithm without the duality gap policy.This table also shows that the average decrease in computational time is now 72.91% from the original CPU times, and it has decreased a further 50.05%from this new computational time when incorporating the duality gap with the enhancements.This is a promising result and demonstrates that the CPLEX adaptive policy has a large and positive effect on the overall efficiency of this enhanced algorithm.
It is important to recognise that for some values of p , such as p = 10 , the total duration could be slightly increased as in this instance the majority of time is spent computing Z -maximal circles.This is because in the first iteration, we do not know whether the majority of time will be spent on computing the Z -maximal circles or solving the problem in CPLEX as CPLEX has not run yet.To respond to this issue, we have therefore set a duality gap of 0.5% for the first iteration.

Overall computational results
Our algorithm was tested on the TSP-Lib data sets rat 575, rat 783, pr 1002 and rl 1323.For information, the data set rat 783 represents a 783-rattled grid problem, and the data sets pr 1002 and rl 1323 refer to a 1002 and 1323-city problem respectively.As we aim to obtain optimal solutions, we used the best known heuristic results from Elshaikh et al. (2015) as our initial upper bound.This deviates from the method previously used, where the initial upper bound was found using the simple H 2 heuristic whose solutions may be relatively loose and hence may require an unnecessarily larger overall computational time.Note also that the computational times given here do not include this heuristic step, but these times are recorded in Elshaikh et al. (2015) .
As these data sets are very large, a maximum time limit of 24 hours was set for each value of p .If the algorithm happens to take longer than the cutoff time, the program is terminated and the upper bound at that time is recorded as the best feasible solution.
Tables 7 -10 are arranged similarly to the tables in Section 5 , with the newly found optimal solutions highlighted in bold.However, extra information for the computational time spent in CPLEX is provided.In order to establish how much computational time cannot be improved on (the last iteration) the column representing the time spent in CPLEX is now divided into two, with one half showing the total time spent in CPLEX and the other half showing how long the last iteration took in CPLEX.Therefore, in the instance where the algorithm reaches the maximum time limit, the result in the second half of this column may not be showing the time spent to reach optimality.However, in each of these circumstances, no further feasible solution was found in the final iteration (except for the case where n = 783 , p = 40 ).Thus, this   indicates that the solution found in the previous iteration may be the optimal solution.Furthermore, in the instance where n = 783 and p = 40 , a feasible solution was found but the duality gap policy value had not been reached.The program was therefore allowed to run for a further hour (with the solution found at this iteration as its new upper bound) to see if this solution could be improved.Again, no further feasible solution was found which shows that the last feasible solution could be optimal.This last feasible solution found is the one given in Table 8 .
It is important to note that for smaller values of p (i.e.p = 10 for pr 1002 and p ≤ 20 for rl 1323) computer memory becomes an issue leading to no results being found.This could be due the initial upper bound being higher in these instances, leading to a relatively large number of circles being considered and thus making the ILP model too big to be handled.
In summary, the results show that the revised Drezner optimal algorithm can now find very good and even optimal solutions for these large data sets.In addition, we can also claim that optimal solutions are found for the first time for the large data sets such  a This excludes computational time for the heuristic step.
* best feasible solution found within 86400 seconds.+ could not be computed due to computer memory.
⊥ no feasible solution found in the last iteration within the time limit allowed.
as n = 575 , n = 1002 and n = 1323 and some for n = 783 while requiring a reasonable amount of computational time only for such strategic decision problems.

Conclusions and suggestions
This paper has revisited an optimal and interesting algorithm proposed by Drezner (1984) thirty years ago to solve the continuous p -centre problem.Opportunities to improve the algorithm were highlighted, and enhancements were developed, mathematically supported and empirically tested.The two areas of interest include the way the Z -maximal circles are identified from one iteration to the next, and the proposed adaptive CPLEX scheme to find a compromise solution at each iteration between the quality of the feasible solution and the optimal solution when solving the covering problem CP 1 .
The proposed algorithm was tested on five existing TSP-Library data sets, namely pr 439, rat 575, rat 783, pr 1002 and rl 1323 for p = 10 , . . ., 100 .The results show that the enhanced optimal method gives a very significant decrease in computational time which sometimes reaches an average reduction of 96%, yielding an algorithm that is superior, faster and more efficient meaning that it can be used to optimally solve the continuous p -centre problem for large data sets for the first time.
One potential research avenue which we believe to be useful would be to incorporate a fast and good heuristic to generate a feasible solution to the covering problem CP 1 instead of using CPLEX all the time.However, as mentioned earlier, at a certain iteration CPLEX or equivalent commercial solver needs to be used to prove infeasibility as this task is mandatory and cannot be performed by a heuristic to guarantee infeasibility.This leads to adopting a new strategy that could combine the exact method and the heuristic approach to solve CP 1 which would identify the appropriate time when the switching from using the heuristic to CPLEX should take place.This is a challenging but interesting task that deserves a thorough investigation.Lastly, research issues related to the tightening of the checking area and in the way the demand points are recorded during the search could also be worth enhancing even further.These aspects are currently being investigated.

Fig. 5 .
Fig. 5. Comparison on CPU time for the enhancements.

Table 1
Number of Z -maximal circles required & previously identified for the first 10 iterations ( n = 439 , p = 100 ).
a This excludes computational time for the H 2 heuristic.

Table 5
CPLEX Durations (seconds) for both the total and the last iteration in the case of n = 575 TSP-Lib

Table 6 n
= 575 TSP-Lib with enhancements and duality gap policy.
a This excludes computational time for the H 2 heuristic.

Table 7
Solutions for n = 575 TSP-Lib using the revised Drezner's algorithm starting from best heuristic value.
a This excludes computational time for the heuristic step.

Table 8
Solutions for n = 783 TSP-Lib using the revised Drezner's algorithm starting from best heuristic value.

Table 9
Solutions for n = 1002 TSP-Lib using the revised Drezner's algorithm starting from best heuristic value.This excludes computational time for the heuristic step.+ could not be computed due to computer memory. a

Table 10
Solutions for n = 1323 TSP-Lib using the revised Drezner's algorithm starting from best heuristic value.