fcmaes - a Python 3 gradient-free optimization library

Traditionally, even when heuristic methods are used, the optimizer is designed for one specific problem. See UAV1, for example, where three heuristic approaches, GA, ACO, and PSO, are compared. We already showed in a fork, UAV2, that results may improve when we give up this idea and apply a generic algorithm instead.

Packing, scheduling, and vehicle routing all have many practically relevant variants. Adapting and reimplementing successful algorithms for each new variant is demanding. In vehicle routing, more generic approaches already exist that cover multiple variants: A generic exact solver for vehicle routing and related problems. Better still would be to avoid creating a problem-specific algorithm entirely.

In that case, we only need an efficient implementation of the objective function, the "scoring" of a solution. That part can often be adapted more easily to many variants. Such a "one for all" approach is useful in several situations:

To study the tradeoff, "how much do we lose with the generic approach", we can test it on well-known problem variants where researchers had decades to improve specialized algorithms. MMKP and VRPTW are good candidates because they are important in practice and heavily studied. Still, we do not see these benchmark problems as the main application area of the approach.

CPU vendors are investing in many-core architectures, and cloud computing offers affordable large-scale parallelism. This trend extends the reach of continuous optimization far beyond its traditional application areas. In this setting, optimization algorithms need different tuning. High diversity across repeated runs becomes much more important. This may explain why the CR-FM-NES → BiteOpt sequence performs better than the alternatives we tried.

Note: CR-FM-NES does not have this property in general. We need to lower the initial step size parameter input_sigma.

Note: The fcmaes implementation of parallelism works much better on Linux. On Windows, consider using the Linux subsystem for Windows WSL.

We have n sets composed of mutually exclusive items. The goal is to select exactly one item from each set and maximize the overall value without violating a family of resource constraints. For each resource, the sum of the resource values of the chosen items must not exceed a resource-specific threshold.

See MMKP-Ants for the mathematical details.

Note: Hard multidimensional multiple choice knapsack problems states: "An instance is hard to solve when all classes contain the same profit vector and the weights are correlated with the profits". The selection of benchmark problems may influence the results.

We use the well-known MMKP benchmark problems "I01-I13", which are cited in most related publications. The literature reports excellent results (see hybrid algorithm), but we did not find the corresponding solution instances. We therefore applied Simulated Annealing for 1E8 iterations with multiple retries to produce at least some kind of "reference solutions", which can be verified:

Note that both hybrid algorithm and Hifi report better solutions. The latter used a 250 Mhz CPU with 128 Mb of RAM and generated an I13 result of 98429 in 160 seconds single-threaded. Please contact me if you know where the corresponding code can be found.

These benchmark results for continuous optimization were computed on a 16 core AMD 5950x CPU using 32 parallel threads with the code at mmkp.py:

As you can see, we lose about 0.5% accuracy compared to the reference solutions above. We still reach 2% accuracy in less than a minute, even for the larger instances. The following diagram shows the relation between the number of groups or decision variables and the time to reach 2% accuracy:

The relation is almost linear. That suggests we can handle even larger instances with reasonable effort.

Github repositories related to MMKP are:

Both hybrid algorithm and Hifi report better results, but there seems to be no related open source code available.

The complete code for the MMKP problem is at mmkp.py. To apply continuous optimization, we represent a problem instance as Python class MMKP:

We extract the problem parameters by parsing the instance file:

The number of decision variables, self.dim, is equal to the number of groups. The boundaries are [0, self.l-1E-12]. The fitness function maps each continuous decision vector to a vector of integers in the [0, self.l-1] interval using numpy’s astype(int). These integers represent the selected item for each group. fitness delegates to a fast numba function, fitness_, which checks resource consumption and returns the overall value together with a penalty value for resource violations.

This function is much easier to adapt to problem variants than optimization algorithms that use the "internal structure" of a problem instance. Those algorithms usually support incremental changes of a solution by computing the score delta. See for instance saMultiChoiceKnapsack.cpp.

As an example, assume we want not only to maximize the sum of item values, but also to achieve balanced resource consumption. We can express this by the standard deviation of the consumed resources divided by their availability. Only a small modification of the fitness or objective function is required:

We applied a specific weight 10 to the standard deviation. For smaller problem instances, we could instead apply a multi-objective algorithm to generate the complete pareto front. That algorithm also comes with fcmaes and is used in many other tutorials.

Apply the same modification to saMultiChoiceKnapsack.cpp. Hint: This may be a bit tricky. Note that this optimization algorithm is not only problem-specific, it is benchmark-specific. It exploits the fact that for most benchmarks the resource limit for all resources is equal. Here double ratio = value/weight the resource consumption weight is not normalized using the resource availability as it should be.

Parallelization of optimization runs and of the optimization algorithm comes for free if we use the fcmaes library:

These settings are sufficient if you aim for a 2% deviation from the optimum (stop_fitness = 2.0). A 1.5% deviation is much harder to achieve, so you may want to increase the number of retries, num_retries. The number of parallel retries depends on the CPU. For the AMD 5950x it is 32, and this can be overwritten using the workers argument. We experimented with smaller population sizes for smaller instances, but there is not much to gain. Both CR-FM-NES and BiteOpt are largely self-adapting.

We plan routes for a fleet of vehicles to serve a set of customer demands. Each vehicle has a capacity constraint on the total demand it serves. Each customer node has a time window, and customer-specific service times must be considered. Possible objectives are the number of vehicles used and the overall distance traveled by all vehicles. If only one vehicle is available, the problem becomes a variant of TSP, the Traveling Salesman Problem.

Possible variants include variable vehicle speed, noisy distances or demands, and additional constraints.

Many different benchmarks are used in the literature. We choose the 100-customer instances of Solomon’s benchmark http://web.cba.neu.edu/~msolomon/problems.htm because reference solutions are available and the benchmark is cited in most related publications.

There are two different objectives for Solomon’s VRPTW benchmark:

The single-objective variant can be solved almost perfectly using or-tools, so we choose that one. See optimize_or.py for the or-tools implementation used to generate the reference results for comparison here. We found other reference solutions at galgos, but some of them did not pass our validation. These solutions assume rounding of the distances, which makes them incompatible with the interpretation of the problem used here.

Besides optimize_or.py, there are many implementations for this problem. http://vrp.galgos.inf.puc-rio.br/index.php/en/links links to some of them. That is not the case for optimization without a problem-specific algorithm. Maybe that is because most continuous optimization algorithms do not work well here.

The complete code for the VRPTW problem is at vrptw.py. To apply continuous optimization, we represent a problem instance as Python class VRPTW:

We extract the problem parameters by parsing the instance file:

The number of decision variables, self.dim, is equal to the number of customer locations. Note that all tables include a 0-entry for the start and end location. We use the [0, 1] interval as boundary. The fitness function converts the continuous argument vector into a list of unique integers using np.argsort(x). fitness then delegates to a fast numba function, fitness_, which executes all tours and evaluates both the objectives and the constraints.

Instead of ignoring the vehicle number, we include it with some weight, even though we are ultimately interested only in the overall distance. This acts as a heuristic that supports the optimization process. We observed that good solutions usually use a fairly low number of vehicles.

You will probably notice that this code is almost exactly the same as for MMKP above. Even the popsize parameter is the same.

Modify the fitness function to handle a problem variant with noisy demands. For example, when using demands[customer], multiply it with a random factor in the [0.8,1.2] interval. You then need to call fitness_ multiple times and use a "worst case" value as the fitness. See link:TSP.html where this was done for TSP. Can you do the same using or-tools by modifying optimize_or.py?

Continuous optimization is performed by a sequence of CR-FM-NES and BiteOpt, with 64 runs in total and 32 runs performed in parallel. On an AMD 5950x 16 core CPU, this takes about 7 minutes.

Compared to the or-tools result, which serves as the reference, we lose about 0.2% on the clustered problem instances and about 3% on the random problem instances.