To refine the facility distribution computed in Phase I, we used a
limited memory BFGS-method (L-BFGS) implemented by Zhu [27].
We also experimented with a specific implementation of a
quasi-Newton method [9] and a truncated Newton method [15],
but our computational experience suggests that the L-BFGS code is faster and
more reliable on our problem than the two other programs.
Since details about the implementation used can be found
elsewhere [27, 4, 5], we discuss here only
changes introduced in order to maximize the performance of the code on our
specific problem.
Here, 
and 
are the function values from the current and the last step,
is the gradient at the current point projected on the
(2N-3)-dimensional cube defined by the constraints and eps is the
machine precision.
To compute a solution of the facility dispersion problem with
highest possible accuracy, we implemented an exact gradient evaluation and used
the following stopping criterion. The code stops when
or when
or when the line search was unsuccessful. (The line search consists of
successive polynomial interpolation and tries to enforce certain
curvature conditions. A line search is deemed to be unsuccessful if there are
more than 10 function evaluations during the search.) We decided to use the first stopping
criterion (17) because of our high accuracy demand. This
criterion is an extremely conservative one in the sense that rounding
errors will usually prevent the inequality from becoming true.
Our computational results suggest that only one run of the simulated annealing algorithm, followed by one run of the L-BFGS method is necessary to produce a good approximation to a global optimum. Moreover, weaker stopping criteria than the ones described above usually led to integration formulae which do not perform as well on test functions as formulae computed with stricter termination rules.