To refine the facility distribution computed in Phase I, we used a limited memory BFGS-method (L-BFGS) implemented by Zhu . We also experimented with a specific implementation of a quasi-Newton method  and a truncated Newton method , 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 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.