fcmaes - a Python 3 gradient-free optimization library

If you just want to experiment with QD algorithms, start with the simple example arm.py, which simulates a planar robot arm. Change the parameters or algorithms to get a feel for how they behave.

In MapElites.adoc we looked at QD optimization from the perspective of an extremely hard problem to understand its limits. We found that QD methods are still applicable even on very rugged solution landscapes such as Cassini2, provided that:

These methods need many fitness evaluations. That raises an obvious question: what happens when each evaluation is expensive because it runs a detailed simulation?

Technical simulations are usually used to understand how input parameters affect system behavior and to identify good parameter settings. In practice, "good" rarely means only one objective. Constraints matter too, and some of them may be flexible. If a small relaxation of a constraint leads to a large improvement, we may decide that the tradeoff is acceptable. For that reason, it is often not a good idea to encode every constraint directly into the optimization problem from the start. It is usually better to explore the design space first and decide later which constraints are acceptable. The same applies to multiple objectives. We can postpone the choice of their weights until we have seen the available tradeoffs.

This means moving away from the usual idea of a single objective, or even a fixed set of objectives, as required by standard optimization algorithms. One option would be to sweep the parameter space with a grid or random samples. That can work for small problems, but it breaks down as the number of decision variables grows. More importantly, the tessellation should usually not be based on the input parameters. It should be based on derived properties that describe the behavior of the simulation. For example:

Imagine a power plant simulation where we care about efficiency, defined as power divided by heat. Then the tessellation should be based on the derived properties power and heat, not on the input parameters. This leads to two questions:

The first question is straightforward. Behavior is observable during the simulation, so the fitness function can return an additional behavior vector. CVT Map Elites proposes a Voronoi tessellation that can be implemented efficiently with sklearn.neighbors.KDTree and sklearn.cluster.KMeans. We tessellate the behavior space into niches and store the elite of each niche, meaning the fittest solution found there, in a QD archive. Building the KDTree can be expensive if we use 10000 or more niches, but caching the niche means helps. If the behavior space is normalized to [0,1] in each dimension, these cached means can even be shared across problems as long as the number of dimensions and niches stays the same. fcmaes does this caching automatically.

The second question is harder. We want to do better than random sampling and focus on interesting parts of the solution space. So we reintroduce an objective, not to find one global optimum, but to define what "interesting" means. There are different ways to handle multiple objectives. For example, we could use a weighted sum as the fitness objective and still preserve diversity by using the individual objectives as behavior dimensions. In the rest of this tutorial we focus on the simpler case of a single fitness objective. In the power plant example, "interesting" means "efficient", so efficiency becomes the fitness objective.

qdo_yahpo provides a framework for testing QD hyperparameter optimization. It uses a surrogate model that can be executed single-threaded, which makes parallel testing easy. We applied these benchmarks to fcmaes diversifier in yahpo.py.

The main differences are:

Being able to choose the improvement-emitter algorithm is important, and the fcmaes diversifier meta-algorithm supports that directly. The problems discussed here work best with either CR-FM-NES or CMA-ES as the emitter. In MapElites.adoc, which applies the fcmaes diversifier to a space mission design problem through elitescass2.py, the best choice is a DE → CMA-ES sequence. In Scheduling.adoc (section QD update) it is BiteOpt. See also https://github.com/google/evojax/pull/52, where I created a PR for EvoJax that introduces the same idea in the machine learning domain.

fcmaes diversifier also performs very well on QD hyperparameter optimization, although a direct comparison is difficult because fcmaes uses Voronoi tessellation while qdo_yahpo uses a grid.

If a surrogate model is available, as in yahpo.py, parallelization is much easier because there is no GPU bottleneck. Without a surrogate, hyperparameter optimization often depends on a GPU or TPU that cannot easily be shared, which usually limits optimization to one thread per device.

The complete code for this example is here: powerplant.py. PowerPlant.adoc explains the single-objective and multi-objective variants. Here we extend the same example with QD methods.

The power plant simulation is based on tespy, a Python framework for thermal engineering systems. We vary the pressure at two extraction connections. These two pressures are the decision variables. After each simulation run, we divide power by heat to compute the efficiency that we want to maximize.

The QD behavior vector desc stores power and heat separately. Note that calculate_qd is protected by with threadpoolctl.threadpool_limits(limits=1, user_api="blas") so the simulation runs single-threaded. That avoids interference with the parallel optimization.

The QD optimization is executed with diversifier.minimize. It runs MAP-Elites ('solver':'elites') and a CMA-ES improvement emitter ('solver':'CMA_CPP') in parallel, with half of the available threads assigned to each. qd_bounds normalizes the behavior values. max_evals=25600 limits the total number of fitness evaluations, while 'max_evals':200 limits one improvement-emitter run. In this example, CMA-ES is the best emitter base algorithm. That is not the case for most of the other simulation-based problems discussed later.

The resulting diagram shows how efficiency is distributed across different power and heat values. It becomes easy to identify the most efficient solutions for specific power and heat limits. Because fcmaes can store and reload the QD archive, you can make these decisions later in a separate analysis step. You can also restart the optimization from a stored archive with different optimization parameters. It is even possible to change the number of niches or the definition of the behavior vector in between, although that requires recomputing fitness values for the stored solutions.

The complete code for this example is here: vilar.py. Sweep.adoc explains the single-objective and multi-objective versions. Here we add QD methods.

In Mechanisms of noise-resistance in genetic oscillators, Jose M. G. Vilar introduced a biochemical model of a circadian clock that lets organisms keep an internal sense of daily time. The model can be simulated with GillesPy2; see VilarOscillator.py. The Vilar model has 15 parameters, and we want to answer these questions:

We use scipy’s argrelextrema to identify the maxima of species R. Then we compute the standard deviation of the amplitudes and of the distances between peaks. Small values for these standard deviations indicate a steady oscillation, so we use them as objectives. ws = sdev_peak_dist/3.0 + sdev_amp/30.0, the normalized weighted sum of these standard deviations, serves as the fitness value. We also include the frequency in the behavior vector to increase diversification.

This time we use CR-FM-NES as the base algorithm for the improvement emitter ('solver':'CRMFNES_CPP') and run it in parallel with MAP-Elites ('solver':'elites'). We choose a small population size because the simulation is expensive. Even with parallelization, we reach only about 8 simulations per second.

The resulting diagram also includes a second run where we maximize the objective ws = 2 - (sdev_peak_dist/3.0 + sdev_amp/30.0). Both diagrams look quite similar, which suggests that the objective plays only a small role here.

The complete code for this example is here: crypto.py. CryptoTrading.adoc explains the single-objective and multi-objective versions. Here we add QD methods.

When we optimize a trading strategy on historical data, the main risk is overfitting to the past. The example already reduces that risk by optimizing ROI (return on investment) across 4 tickers and using the geometric mean ROI as fitness. This ROI is normalized against the HODL ROI, meaning the return we would get by buying and holding the entire time. We can now use the 4 normalized ROI factors as the behavior vector and generate a diverse set of solutions.

Once we have this archive of diverse solutions, we can count how often each parameter value occurs.

This gives us a more reliable indicator of which parameter values work well:

The complete code for this example is here: lrgv.py. See also Water.adoc.

The Lower Rio Grande Valley (LRGV) framework models a risk-based water supply portfolio problem. A single city must choose an efficient mix of market-based sources and traditional reservoir sources while minimizing the risk of not having enough water available when needed. An option-based market lets the city secure water later at a fixed price by paying an option price in advance.

We forked the original code at https://github.com/dietmarwo/LRGV so it can be called from Python and made it reentrant. That greatly increases the number of simulations we can execute per second, so running 400000 fitness evaluations becomes easy.

We configure the framework with five objectives:

For QD, we use a weighted sum of these five objectives as the fitness value. The behavior vector contains the five objectives themselves.

Since there are five objectives, each diagram shows three of them.

There are not many complex real-world multi-objective benchmark problems in the public domain. One of them is the Mazda Benchmark Problem, developed jointly by Mazda Motor Corporation, the Japan Aerospace Exploration Agency, and Tokyo University of Science. It has 222 discrete decision variables, 54 inequality constraints, and multiple objectives. The task is to design three cars at the same time while minimizing weight and maximizing the number of shared sheet-thickness parts, which reduces production cost. The original constraints represent expensive collision simulations used to assess car safety. In the benchmark, these simulations are replaced by response-surface approximations, which act as a domain-specific surrogate model. After generating solutions on the approximate model, the real collision simulations can be applied later to filter out solutions that do not hold up in the real world.

In Surrogate we described how to solve this problem with MODE, the fcmaes multi-objective algorithm.

With QD, the Mazda benchmark now yields more than 7000 diverse solutions, far more than a typical MO algorithm returns.

The figure above shows progress over time:

The constraints approximate expensive physical simulations, namely crash tests. Satisfying the approximate constraints does not guarantee that the real constraints are also satisfied. The opposite can also happen: the approximation may be more restrictive than reality.

So not all 7000 solutions will be equally interesting. The most valuable ones may be the few hundred solutions near the border between valid and invalid regions. These can be the starting point for more detailed follow-up studies with more expensive simulations that verify the constraints.

There is still little literature on applying QD methods to expensive simulations, so there is not much we can compare against. These results are therefore preliminary, but they can still serve as a useful baseline for future work.

We found that: