fcmaes - a Python 3 gradient-free optimization library

Solving the ESP2, ESP3 and PitzDaily Fluid Dynamics Benchmark Problems

This tutorial

Shows how to solve the ESP and PitzDaily real-world problems using fcmaes differential evolution.
Compares the results with the surrogate-model-based solver MVRSM.

Motivation

From Black-box Combinatorial Optimization using Models with Integer-valued Minima:

"When a black-box optimization objective can only be evaluated with costly or noisy measurements, most standard optimization algorithms are unsuited to find the optimal solution. Specialized algorithms that deal with exactly this situation make use of surrogate models."

ESP2 is exactly this kind of problem. The practical question is whether surrogate-model-based solvers really have a clear advantage here.

Evaluating the objective function does three expensive things:

It loads a docker image containing the open source fluid dynamics software OpenFOAM.
It simulates an electrostatic precipitator configured by the objective vector.
ESP2 uses both 49 continuous and 49 discrete variables.

An electrostatic precipitator is a large gas filtering installation. Its efficiency depends on how evenly the incoming gas is distributed. The installation contains slots of different types, and each slot type affects the flow distribution differently. The goal is to find the configuration that produces the best distribution.

How to apply fcmaes-DE

fcmaes-DE has two properties that make it attractive for very expensive optimization problems:

It converges quickly because it is based on the DE/best/1 strategy.
It supports parallel function evaluation.

The first experiments with fcmaes-DE on ESP2 and ESP were disappointing. For scheduling and routing problems such as job shop, GTOC11, or noisy TSP, we had successfully used the np.argmin or sorting trick to map continuous variables to unique integer values. ESP2 is different. Its discrete values do not need to be unique, and their range is very small: [0,3].

At that point it became clear that fcmaes needed direct mixed-integer support. That support was then added to fcmaes-DE and fcmaes-MODE, the new multi-objective and constraint algorithm, in both the Python and C++ implementations.

You can now tell the algorithm which variables are discrete by using the ints parameter. For example, ints = [True, True, False] means that the first two variables are integer variables.

from scipy.optimize import Bounds
from fcmaes import decpp, de

problem = ESP2

ret = decpp.minimize(problem.evaluate,
            dim = problem.dims(),
            bounds = Bounds(problem.lbs(),problem.ubs()),
            popsize = 31,
            max_evaluations = 10000,
            workers = 12,
            ints = [v != 'cont' for v in problem.vartype()]
            )

print(ret.x, ret.fun)

For logging, you can wrap problem.evaluate. Because the objective is evaluated in parallel, use mp.RawValue to monitor progress:

        self.evals = mp.RawValue(ct.c_long, 0)  # writable across python processes
        self.best_y = mp.RawValue(ct.c_double, math.inf) # writable across python processes

On Windows, replace decpp.minimize with de.minimize, the Python variant. Parallel retry works on Windows, but parallel function evaluation has issues.

On an AMD 5950x 16-core processor, we get about 0.7 ESP2 evaluations per second. Docker calls do not scale perfectly, but the speedup is still useful.

Lowering popsize to 16 or 8 speeds up convergence, but it also greatly increases the chance of getting stuck in a local minimum.

The 8 MVRSM runs were performed in parallel on a 16-core machine. The DE runs use parallel function evaluation instead. MVRSM either needs some luck, or you can interpret the 8 parallel runs as a crude form of parallelization where only the best result matters. For a fair comparison, only the best MVRSM run should be considered.

There is not enough data here to declare a winner. Both approaches reach excellent results, unlike all other optimization methods tested, including fcmaes-DE without mixed-integer support.

Still, these results make it hard to justify a general claim that surrogate-based methods are superior.

How is mixed integer support implemented

The mixed-integer support is not complicated. It is a small but effective modification.

The first change affects the bounds. Internally, the algorithm still works with continuous variables. For each discrete choice, we therefore assign the same interval width inside the bounds. In practice, we add 0.99999 to the upper integer bound and assume that the fitness function truncates to an integer.

The second change affects mutation and crossover. Generated values are truncated to the next integer value.

This matters because candidates are created from differences between vectors. That is why the method is called Differential Evolution. We want those differences to stay discrete for integer variables.

We also add random mutations on the integer variables. The mutation probability depends on the number of integer variables.

This was motivated by observed stagnation in the integer variables during later optimization stages. Looking at the MVRSM code, I noticed that it applies many mutations to discrete variables, so trying something similar seemed reasonable. If you use the Python variants, you can define your own sample modifier and override the default behavior. Please let me know if you find an improvement.

    # default modifier for integer variables
    def _modifier(self, x):
        x_ints = x[self.ints]
        n_ints = len(self.ints)
        lb = self.lower[self.ints]
        ub = self.upper[self.ints]
        min_mutate = 0.1
        max_mutate = 0.5 # max(1.0, n_ints/20.0)
        to_mutate = self.rg.uniform(min_mutate, max_mutate)
        # mututate some integer variables
        x[self.ints] = np.array([x if self.rg.random() > to_mutate/n_ints else
                           int(self.rg.uniform(lb[i], ub[i]))
                           for i, x in enumerate(x_ints)])
        return x

ESP3

ESP3 is another mixed-integer variant of the ESP CFD problem, but with fewer continuous variables.

from scipy.optimize import Bounds
from fcmaes import decpp, de

problem = ESP3

ret = decpp.minimize(problem.evaluate,
            dim = problem.dims(),
            bounds = Bounds(problem.lbs(),problem.ubs()),
            popsize = 24,
            max_evaluations = 5000,
            workers = 12,
            ints = [v != 'cont' for v in problem.vartype()]
            )

print(ret.x, ret.fun)

For this problem, popsize can be reduced to 24. Parallel execution on an AMD 5950x CPU gives an execution time of about 1.5 sec per evaluation.

In Bliek21 you can find results for surrogate-based optimizers on this problem.

PitzDaily

PitzDaily is another benchmark included in ESP and PitzDaily. It studies how combustion affects mean flowfield properties such as mixing layer growth, entrainment rate, and reattachment length. PitzDaily provides a useful visualization of the problem. It was chosen as an OpenFOAM tutorial because the setup is relatively simple: only continuous variables and a low dimension of 10. fcmaes-DE can solve it easily.

from scipy.optimize import Bounds
from fcmaes import decpp, de

problem = PitzDaily

ret = decpp.minimize(problem.evaluate,
            dim = problem.dims(),
            bounds = Bounds(problem.lbs(),problem.ubs()),
            popsize = 24,
            max_evaluations = 5000,
            workers = 12,
            )

print(ret.x, ret.fun)

Here, popsize can also be reduced to 24. No ints parameter is needed because all variables are continuous.

After about 600 seconds, 12 out of 13 runs reach 0.08. The last run also succeeds after about 1100 seconds. One surprising detail is that solutions < 0.079 seem hard to find in the literature. Parallel execution on an AMD 5950x CPU gives about 0.7 sec per evaluation.

ESP4 comparison with Bayesian Optimization and MVRSM

From Frazier2018: "Bayesian optimization is an approach to optimizing objective functions that take a long time (minutes or hours) to evaluate. It is best-suited for optimization over continuous domains of less than 20 dimensions, and tolerates stochastic noise in function evaluations. It builds a surrogate for the objective and quantifies the uncertainty in that surrogate using a Bayesian machine learning technique, Gaussian process regression, and then uses an acquisition function defined from this surrogate to decide where to sample."

We use ESP4 to test whether Bayesian Optimization also applies at higher dimension. ESP4 has dimension 54. Bayesian Optimization is interesting here because:

Some implementations support mixed-integer problems, such as https://github.com/SheffieldML/GPyOpt and https://github.com/wangronin/Bayesian-Optimization.
Some implementations support parallel function evaluation, such as those two and https://github.com/wujian16/Cornell-MOE .

Unfortunately, https://github.com/wangronin/Bayesian-Optimization failed on the docker-based ESP4 problem with:

_pickle.PicklingError: Could not pickle the task to send it to the workers.

https://github.com/SheffieldML/GPyOpt works, but after the initial population it falls back to a single thread. https://github.com/wujian16/Cornell-MOE does not support mixed-integer problems.

So for this comparison we used https://github.com/SheffieldML/GPyOpt in single-threaded mode. To make use of the CPU, we ran 8 optimizations in parallel.

GPyOpt needs about 17.6 sec per evaluation (8 parallel optimizations), so 1000 evaluations need about 17500 sec.
MVRSM needs about 17 sec per evaluation (8 parallel optimizations), so 1000 evaluations need about 17000 sec.
fcmaes DE needs about 1.5 sec per evaluation (12 parallel function evaluations), so 1000 evaluations need about 1500 sec.

Here are the results for fcmaes DE using the mixed-integer enhancement, limited to 5000 evaluations:

Even if we compare DE after 1000 evaluations / 1500 sec with BO after 1000 evaluations / 17500 sec, fcmaes DE with mixed-integer support is clearly ahead. That conclusion does not depend on the higher algorithmic overhead of Bayesian Optimization or on the fact that we could not get parallel evaluation working. Even if we compare only by evaluation count, Bayesian Optimization does not work well for ESP4 at dimension 54.

MVRSM is more competitive against DE when you compare the result after 1000 evaluations / 17000 sec. It is also worth noting that after only 200 evaluations, BO and MVRSM are very close, and both beat Differential Evolution. For problems so expensive that only 200 evaluations are feasible, BO and MVRSM could be the better choice, provided parallel evaluation is supported and actually works. After 1000 evaluations, only DE and MVRSM can be recommended for high-dimensional problems.

Conclusion

To summarize:

We could confirm the superiority of surrogate-based methods like BO and MVRSM for complex mixed-integer CFD simulation problems if only about 200 function evaluations are affordable. But the results after only 200 evaluations are still poor, and parallel execution can increase the feasible evaluation budget.
fcmaes Differential Evolution, thanks to its mixed-integer support, is a serious competitor if more than 500 evaluations can be performed. It also already supports parallel function evaluation.
MVRSM is, like DE, superior to Bayesian Optimization for higher evaluation budgets.
The fcmaes multi-objective solver (MO-DE) with mixed-integer support is ready to be tested in this application area.