The ability to perform per-epoch sensitivity analysis and dispersion index computation means that the EO algorithm parameters will differ from to epoch. We need some general way to record the parameters of the EO algorithm used for each epoch.

EHVI = Expected Hypervolume Improvement

References:
Efficient Computation of Expected Hypervolume Improvement Using Box Decomposition Algorithms. https://arxiv.org/pdf/1904.12672.pdf
Efficient computation of the search region in multi-objective optimization. EJOR 2017.
A Box Decomposition Algorithm to Compute the Hypervolume Indicator. https://arxiv.org/abs/1510.01963

https://arxiv.org/abs/1612.01474

The tournament selection routine creates a vector of probabilities for each candidate solution, and then uses random.choice to select candidates for the MOEA pool selection. However, it does not ensure the probabilities sum to 1, which occasionally leads to errors.

dmosopt needs a more general approach to providing optimizer parameters.

The refactoring to generator-based code has broken the hot start logic. A corresponding unit test should be added as well.

This feature will add a simple multi-objective trust region optimizer, which can be used as a local optimizer after an epoch with a global optimizer.

After the recent refactorings, constraint handling is not applied after initial sampling epoch, but is applied when loading saved evaluations.

Example debug.py script to reproduce:

import logging
import numpy as np
from dmosopt import dmosopt

logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)


def zdt1(x):
    num_variables = len(x)
    f = np.zeros(2)
    f[0] = x[0]
    g = 1.0 + 9.0 / float(num_variables - 1) * np.sum(x[1:])
    h = 1.0 - np.sqrt(f[0] / g)
    f[1] = g * h
    return f


def obj_fun(pp):
    param_values = np.asarray([pp[k] for k in sorted(pp)])
    res = zdt1(param_values)
    logger.info(f"Iter: \t pp:{pp}, result:{res}")
    return res


def zdt1_pareto(n_points=100):
    f = np.zeros([n_points, 2])
    f[:, 0] = np.linspace(0, 1, n_points)
    f[:, 1] = 1.0 - np.sqrt(f[:, 0])
    return f


if __name__ == "__main__":
    space = {}
    for i in range(30):
        space["x%d" % (i + 1)] = [0.0, 1.0]
    problem_parameters = {}
    objective_names = ["y1", "y2"]
    dmosopt_params = {
        "opt_id": "dmosopt_zdt1",
        "obj_fun_name": "debug.obj_fun",
        "problem_parameters": problem_parameters,
        "space": space,
        "objective_names": objective_names,
        "population_size": 50,
        "num_generations": 25,
        "initial_maxiter": 10,
        "optimizer_name": "age",
        "termination_conditions": True,
        "n_initial": 3,
        "n_epochs": 4,
        "surrogate_method_name": None,
        "save": True,
        "file_path": "zdt1.h5",
    }

    dmosopt.run(dmosopt_params, verbose=True)

Traceback (most recent call last):
  File "./MultiObjectiveSurrogateOptimization/debug.py", line 64, in <module>
    dmosopt.run(dmosopt_params, verbose=True)
  File "./dmosopt/dmosopt/dmosopt.py", line 1975, in run
    distwq.run(
  File "/scratch1/08818/fg14/venvs/multiobjectivesurrogateoptimization-py3.9-scratch/lib/python3.9/site-packages/distwq.py", line 2071, in run
    fun(controller, *args)
  File "./dmosopt/dmosopt/dmosopt.py", line 1944, in dopt_ctrl
    dopt.run_epoch()
  File "./dmosopt/dmosopt/dmosopt.py", line 1137, in run_epoch
    strategy_state, strategy_value, completed_evals = self.optimizer_dict[
  File "./dmosopt/dmosopt/dmosopt.py", line 357, in update_epoch
    x_resample = result_dict["x_resample"]
KeyError: 'x_resample'

https://mlopez-ibanez.github.io/irace/index.html

A more accurate implementation of the time limit functionality should include 1) GP model train time measurements and estimates 2) elapsed time and time limit in MOEA optimizer..

The recommended approach for reproducible RNG in recent versions of numpy is to pickle the generator.

This feature will add support for using different optimization methods in successive epochs, i.e. NSGA in epoch 1, followed by a local optimization method in epoch 2, then followed by NSGA again and so on. It will require modifying the various optimizer_method arguments to take an iterable of optimizer method names.

Currently, the predict methods in the Gaussian Process interface module only return the mean at the requested points. However, the variance is also necessary to estimate the uncertainty at those points.

https://github.com/ziatdinovmax/gpax

Currently, the user can specify the objective and initialization methods by simply providing a Python module and object name. However, the optimizer and sampling strategies are hardcoded.

I wonder if it would make sense to change the semantics of the optimizer and sampling_method setting to refer to a Python object instead. For example, "optimizer": "nsga2" would become "optimizer": "dmosopt.NSGA2:NSGA2" etc. In this example, I've used the entry point notation, but for consistency, it might make more to add a optimizer_module option instead.

This would allow users to extend and develop custom strategies by simply pointing to the module. To maintain backward compatibility we could simply keep the current behavior but fall back on importing if the user provides an 'unknown' option.