fcmaes - a Python 3 gradient-free optimization library

Assign service facilities to demand points

This tutorial

shows how to assign service facilities to demand points while accounting for failure of a given number of service points.
compares our approach with https://github.com/netotz/alpha-neighbor-p-center-problem
shows how to adapt the code to a continuous variant of the problem.
See also 5G.adoc, which covers a similar problem.

α-neighbor p-center optimization problem

This tutorial is inspired by https://github.com/netotz/alpha-neighbor-p-center-problem, a new implementation of the discrete variant of the α-neighbor p-center optimization problem. We need to assign multiple service facilities to demand points so the assignment still works when facilities fail.

The goal is to minimize the maximum service distance for all demand points. But we do not use the nearest facility. We use the α-nearest facility instead. This models a worst-case scenario where the α-1 nearest facilities have already failed. There are two variants of this problem:

We choose from a given set of facility locations. This variant has discrete decision variables.
We compute optimal coordinates for the facilities. This variant has continuous decision variables.

Our approach covers both variants.

We will use the BiteOpt optimization algorithm with parallel retry to make good use of modern many-core CPUs.

Motivation

The α-neighbor p-center optimization problem is one of many problems discussed in our tutorials. In several of these domains, specialized algorithms were developed to improve performance.

These algorithms usually exploit properties of the problem domain to update the fitness efficiently for similar solutions, for example by adding, removing, or changing a single service facility. That means they improve a solution through incremental changes only. This has several disadvantages:

Optimization can get stuck more easily in local minima when solutions change only locally and not globally.
Variants of the problem are harder to handle, for example when adding noise, extra constraints, extra objectives, or an objective that also considers the nearest facilities.
The continuous variant, where we search directly for optimal facility coordinates, cannot be handled efficiently by the incremental-change approach.

For small and easy instances, the incremental approach will of course be faster. But a more generic method can recover part of that performance loss by using the parallelism of modern many-core CPUs. As long as the runtime stays feasible, in seconds or minutes, flexibility, ease of implementation, and solution quality often matter more.

Implementation

The complete code for the α-neighbor p-center problem is at anpcp.py. We start with the discrete variant.

We use numba to speed up evaluation of a facility selection stored as the index array selection. The distance matrix between demand points (users) and facilities (distance) is a two-dimensional array, and alpha specifies which nearest facility we use. To avoid sorting all distances, we use np.partition to separate the α-nearest distances in linear time. An efficient fitness implementation is essential because, unlike https://github.com/netotz/alpha-neighbor-p-center-problem, we do not rely on partial fitness updates.

    @njit(fastmath=True)
    def fitness_(selection, distance, alpha):
        selected = distance[:,selection]
        partitioned = np.partition(selected, alpha)
        return max([max(d[:alpha]) for d in partitioned])

    ...
    class ANPCP():
        ...
        def init_json(self, json_file):
            with open(json_file) as json_file:
            ...
                self.fnum = len(self.facilities)
                self.bounds = Bounds([0]*self.dim, [self.fnum-1E-9]*self.dim)


        def fitness(self, x):
            selection = selection_(x.astype(int), self.fnum)
            return fitness_(selection, self.distance, self.alpha)

A problem instance is represented by the Python fitness class ANPCP.

To use a continuous optimizer such as BiteOpt, we need to do two things:

Define bounds for the continuous decision variables. We use the interval [0, fnum[, where fnum, the total number of facilities, defines the number of choices.
Map those variables efficiently to an index array representing the selected facilities. x.astype(int) converts the decision vector x to an integer array, but it may contain duplicates. The numba function selection_(s, n) below fixes that by returning unique indices in the same range.

    @njit(fastmath=True)
    def next_free_(used, p):
        while used[p]:
            p = (p + 1) % used.size
        used[p] = True
        return p

    @njit(fastmath=True)
    def selection_(s, n):
        disjoined_s = np.zeros(s.size, dtype=numba.int32)
        used = np.zeros(n, dtype=numba.boolean)
        for i in range(s.size):
            disjoined_s[i] = next_free_(used, s[i])
        return disjoined_s

Loops like these should not be used in plain Python without numba. We run the optimization as a parallel retry. retry.minimize uses all available CPU cores via mp.cpu_count(), but the number of workers is configurable.

    def optimize(anpcp, opt, num_retries = 32):
        ret = retry.minimize(wrapper(anpcp.fitness),
                                   anpcp.bounds, num_retries = num_retries,
                                   optimizer=opt, logger=logger())

Finally, we choose a concrete problem instance and an optimizer. For this task, BiteOpt combined with parallel retry performs best. For comparison, we use an adapted TSP instance from https://github.com/netotz/alpha-neighbor-p-center-problem.

    anpcp = ANPCP(12, 2) # p = 12, alpha = 2
    anpcp.init_tsp('data/rl1323_993_330_4.anpcp.tsp')
    popsize = 7 + 12*anpcp.dim
    max_evaluations = 300000
    opt = Bite_cpp(max_evaluations, popsize=popsize, M=8)
    optimize(anpcp, opt, num_retries = 32)

Note that BiteOpt is self-adapting, so it also works well without the popsize setting.

Results

When we execute anpcp.py, we get a result in less than one minute on an AMD 5950x 16-core CPU for the rl1323_993_330_4.anpcp.tsp instance with 330 facilities and 993 demand points, selecting 12 facilities with alpha = 2.

36.28 5597295 154280.0 4190.0 [251.26756648242048, 220.01575093780303, ..]
54.94 5460 32 300000 4480.000000 0.00 0.00 [] [101.89788111176522, 329.8952205906099, ...
57.92 165745 32 9600000 4190.000000 4319.69 94.81 [4190.0, 4190.0, 4196.0, 4201.0, ...] [296.44872332608435, 7.36882765692593, ...]
selection =  [296   7  88 162 272  81 133 252 221  53 251 115]
value =  4190.0

The resulting selection [296, 7, 88, 162, 272, 81, 133, 252, 221, 53, 251, 115] has value 4190.0. Repeated runs produce similar results.

Excercise

Compare the performance of different optimization algorithms from fcmaes.optimize such as de_cma, Cma_cpp, De_cpp, Da_cpp, Csma_cpp, Bite_cpp and Crfmnes_cpp.

Comparison

What happens if we try the same problem using https://github.com/dietmarwo/fast-cma-es/blob/master/examples/anpcp/ ?

from models.instance import Instance
from models.solver import Solver

filepath = os.path.abspath("../data/rl1323_993_330_4.anpcp.tsp")
instance = Instance.read_tsp(filepath)
solver = Solver(instance, 12, 2, True)
solver.grasp(30000)

We configure a runtime of 30000 seconds and start 16 runs in parallel, one for each core of the 16-core CPU, because no out-of-the-box parallelization is provided. Even with that large amount of CPU time, the best result we obtained was value = 4388.

We also observed instances where https://github.com/dietmarwo/fast-cma-es/blob/master/examples/anpcp/ was superior, for example large random instances with at least 2000 facilities and users. The open question is how relevant such random instances are for real-world applications.

Locate Service Facilities

There is also a continuous variant of the problem. Here we are not given a fixed set of candidate facilities. We know only p, the number of facilities to place, and we search for optimal coordinates directly. Once we know which regions look promising, we can look for concrete service locations there. That leads back to the first problem variant.

The code is at anpcp.py.

Since we are now using continuous optimization, only small code changes are needed. optimal_algorithms_for_anpcp describes a problem-specific algorithm that gives slightly better solutions, below 1%, for very large instances. That is the price of taking the generic route. But as soon as we add constraints, extra objectives, or noise, adapting the specialized algorithm becomes much harder.

Implementation

The complete code for the continuous variant is at anpcpc.py.

Only minor changes are needed in the objective function. Instead of forwarding a facility selection, we now pass the x- and y-coordinates of the facilities to the fitness function.

    @njit(fastmath=True)
    def fitness_(facilities_x, facilities_y, users, alpha):
        distance = calc_distance_(users, facilities_x, facilities_y)
        partitioned = np.partition(distance, alpha)
        return max([max(d[:alpha]) for d in partitioned])
    ...
    class ANPCPC():
    ...
        def fitness(self, x):
            facilities_x = x[:self.p]
            facilities_y = x[self.p:]
            return fitness_(facilities_x, facilities_y, self.users, self.alpha)

The input vector is split into two halves. One half contains the x-coordinates, the other the y-coordinates.

Results

When we execute anpcpc.py, we get a result in less than 30 seconds on an AMD 5950x 16-core CPU for the rl1323_993_330_4.anpcp.tsp instance with 330 facilities and 993 demand points, selecting 12 facilities with alpha = 2.

27.6 115942 32 3200000 14403864.206926 14766331.71 476172.12 [14403864.21, 14545573.38, 14545573.39, 14545573.39, ...]
facility locations =  [[ 4637.70618771  3245.83435739]
 [ 4547.12658139  3329.12972472]
 [ 9724.74192183  8792.00329984]
 [15344.31817208  2923.37437901]
 [ 9963.97479073  2390.59362575]
 [ 9724.7415082   8792.00297781]
 [15775.90549813  8970.55944954]
 [15237.52041848  2808.93959076]
 [ 3208.16211282  9342.33672938]
 [ 9942.10780989  2659.6540616 ]
 [15918.37951594  8786.078065  ]
 [ 3257.31367395  9307.58761972]]
value =  3795.242312017297

The resulting coordinates have a value of 3795.2. As expected, this is better than 4190.0, the value obtained when choosing from a fixed set of facilities.

Excercise

Again compare the performance of different optimization algorithms from fcmaes.optimize such as de_cma, Cma_cpp, De_cpp, Da_cpp, Csma_cpp, Bite_cpp and Crfmnes_cpp. Note that the results differ significantly from the discrete variant.

Comparison

Compared with the results from α-neighbor p-center optimization problem, we see:

Almost equal results for small and moderate problem instances.
Almost equal results for small facility numbers
Slightly inferior results for large problem instances and facility numbers.

For example, for the pr439_220_219_0.anpcp.tsp instance with 439 facilities, selecting 70 facilities with alpha = 2, we get:

...
600.5 58174583 96877.0 406786.43171511905 [10817.444908968146, 10733.530275534793, ....]
...

This means that after 600 seconds we reach value² = 406786, which gives value = 637.8 after taking the square root. The algorithm optimizes squared distances to save time, so we need to convert the result back. After 600 seconds, no further improvement occurs.

α-neighbor p-center optimization problem reports value = 621.74 after 1888 seconds, which is 2.5% better than our result.

So there is a cost to using a generic algorithm for this variant. Here the specialized algorithm is better. But for most problem instances, the difference is small.

Conclusion

The α-neighbor p-center optimization problem can be solved efficiently by dedicated algorithms.
For most problem instances, the generic approach based on continuous optimization produces equal or better results.
The generic approach often needs more compute resources, but parallelization and an efficient fitness implementation can offset part of that cost.
The BiteOpt algorithm together with parallel retry is a good choice for this problem.