DIFF_EVOLUTION

Overview

The DIFF_EVOLUTION function finds the global minimum of a multivariate function using differential evolution (DE), a stochastic population-based optimization algorithm. Unlike gradient-based methods, DE does not require derivatives and can effectively search large, complex parameter spaces with multiple local minima, making it well-suited for challenging global optimization problems.

Differential evolution was introduced by Storn and Price in 1997 and has since become one of the most widely used evolutionary algorithms for continuous optimization. The algorithm maintains a population of candidate solutions and iteratively improves them through mutation, crossover, and selection operations. For a comprehensive overview, see the Wikipedia article on differential evolution and the SciPy documentation.

The core mutation operation creates trial vectors by combining existing population members. The default best1bin strategy uses the formula:

b' = x_0 + F \cdot (x_{r0} - x_{r1})

where x_0 is the current best solution, x_{r0} and x_{r1} are randomly selected population members, and F is the mutation constant controlling the differential weight. The algorithm then applies binomial crossover with probability CR (the recombination constant) to create the final trial vector.

This implementation leverages SciPy’s differential_evolution optimizer, which supports multiple mutation strategies including rand1bin, best2bin, randtobest1bin, and currenttobest1bin. Key parameters include popsize (population size multiplier), mutation (differential weight, optionally with dithering for adaptive mutation), and recombination (crossover probability). The optional polish parameter enables local refinement using L-BFGS-B after the evolutionary search completes.

For improved global search, use higher popsize and mutation values with lower recombination, though this increases computation time. The algorithm converges when the population’s standard deviation of fitness values falls below the specified tolerance.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=DIFF_EVOLUTION(func_expr, bounds, strategy, maxiter, popsize, tol, mutation, recombination, seed, polish)

func_expr (str, required): Objective expression using x[i] for variable access (e.g., x[0]2 + x[1]2).
bounds (list[list], required): 2D list of [min, max] pairs defining the domain for each variable.
strategy (str, optional, default: “best1bin”): Differential evolution strategy identifier.
maxiter (int, optional, default: 1000): Maximum number of generations.
popsize (int, optional, default: 15): Population size multiplier (population = popsize × dimension).
tol (float, optional, default: 0.01): Relative convergence tolerance.
mutation (float, optional, default: 0.5): Mutation constant in (0, 2]. Pass a 2D list [[min, max]] for adaptive mutation.
recombination (float, optional, default: 0.7): Crossover probability in (0, 1].
seed (int, optional, default: null): Optional integer seed for reproducibility.
polish (bool, optional, default: true): Whether to perform local search polishing on the best solution.

Returns (list[list]): 2D list [[x1, x2, …, objective]], or error message string.

Examples

Example 1: Sphere function with all parameters

Inputs:

func_expr	bounds		strategy	maxiter	popsize	tol	mutation	recombination	seed	polish
x[0]2 + x[1]2	-5	5	best1bin	500	20	0.0001	0.6	0.8	123	true
	-5	5

Excel formula:

=DIFF_EVOLUTION("x[0]**2 + x[1]**2", {-5,5;-5,5}, "best1bin", 500, 20, 0.0001, 0.6, 0.8, 123, TRUE)

Expected output:

Result
0	0	0

Example 2: Shifted quadratic with adaptive mutation

Inputs:

func_expr	bounds		seed	strategy	mutation		recombination
(x[0]-2)2 + (x[1]+1)2	-6	6	321	best2bin	0.5	1.2	0.9
	-6	6

Excel formula:

=DIFF_EVOLUTION("(x[0]-2)**2 + (x[1]+1)**2", {-6,6;-6,6}, 321, "best2bin", {0.5,1.2}, 0.9)

Expected output:

Result
2	-1	0

Example 3: Absolute value without polishing

Inputs:

func_expr	bounds		seed	polish	tol
abs(x[0]) + abs(x[1])	-3	3	999	false	0.000001
	-3	3

Excel formula:

=DIFF_EVOLUTION("abs(x[0]) + abs(x[1])", {-3,3;-3,3}, 999, FALSE, 0.000001)

Expected output:

Result
0	0	0

Example 4: Three-variable optimization

Inputs:

func_expr	bounds		seed	maxiter	popsize
(x[0]-1)2 + (x[1]-0.5)2 + (x[2]+2)**2	-5	5	2024	200	10
	-5	5
	-5	5

Excel formula:

=DIFF_EVOLUTION("(x[0]-1)**2 + (x[1]-0.5)**2 + (x[2]+2)**2", {-5,5;-5,5;-5,5}, 2024, 200, 10)

Expected output:

Result
1	0.5	-2	0

Python Code

import math
import re

import numpy as np
from scipy.optimize import differential_evolution as scipy_differential_evolution

def diff_evolution(func_expr, bounds, strategy='best1bin', maxiter=1000, popsize=15, tol=0.01, mutation=0.5, recombination=0.7, seed=None, polish=True):
    """
    Minimize a multivariate function using differential evolution.

    See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.differential_evolution.html

    This example function is provided as-is without any representation of accuracy.

    Args:
        func_expr (str): Objective expression using x[i] for variable access (e.g., x[0]**2 + x[1]**2).
        bounds (list[list]): 2D list of [min, max] pairs defining the domain for each variable.
        strategy (str, optional): Differential evolution strategy identifier. Valid options: best1bin, best1exp, best2bin, best2exp, rand1bin, rand1exp, rand1binwithmem, rand1expwithmem, rand2bin, rand2exp, randtobest1bin, randtobest1exp, currenttobest1exp. Default is 'best1bin'.
        maxiter (int, optional): Maximum number of generations. Default is 1000.
        popsize (int, optional): Population size multiplier (population = popsize × dimension). Default is 15.
        tol (float, optional): Relative convergence tolerance. Default is 0.01.
        mutation (float, optional): Mutation constant in (0, 2]. Pass a 2D list [[min, max]] for adaptive mutation. Default is 0.5.
        recombination (float, optional): Crossover probability in (0, 1]. Default is 0.7.
        seed (int, optional): Optional integer seed for reproducibility. Default is None.
        polish (bool, optional): Whether to perform local search polishing on the best solution. Default is True.

    Returns:
        list[list]: 2D list [[x1, x2, ..., objective]], or error message string.
    """
    def to2d(x):
        return [[x]] if not isinstance(x, list) else x

    if not isinstance(func_expr, str):
        return "Invalid input: func_expr must be a string."
    if func_expr.strip() == "":
        return "Invalid input: func_expr must be a non-empty string."

    func_expr = re.sub(r'\^', '**', func_expr)

    if not re.search(r'\bx\b', func_expr):
        return "Invalid input: func_expr must reference variable x (e.g., x[0])."

    bounds = to2d(bounds)

    if not isinstance(bounds, list) or len(bounds) == 0:
        return "Invalid input: bounds must be a 2D list of [min, max] pairs."

    processed_bounds = []
    for idx, bound in enumerate(bounds):
        if not isinstance(bound, list) or len(bound) != 2:
            return "Invalid input: each bound must be a [min, max] pair."
        try:
            lower = float(bound[0])
            upper = float(bound[1])
        except (TypeError, ValueError):
            return "Invalid input: bounds must contain numeric values."
        if lower > upper:
            return f"Invalid input: lower bound must not exceed upper bound for variable index {idx}."
        processed_bounds.append((lower, upper))

    dimension = len(processed_bounds)

    valid_strategies = {
        'best1bin', 'best1exp', 'rand1exp', 'randtobest1exp', 'currenttobest1exp',
        'best2exp', 'rand2exp', 'randtobest1bin', 'best2bin', 'rand2bin',
        'rand1bin', 'rand1binwithmem', 'rand1expwithmem'
    }
    if strategy not in valid_strategies:
        return f"Invalid input: strategy must be one of: {', '.join(sorted(valid_strategies))}"

    try:
        maxiter = int(maxiter)
        popsize = int(popsize)
        tol = float(tol)
        recombination = float(recombination)
    except (TypeError, ValueError):
        return "Invalid input: maxiter, popsize must be integers; tol, recombination must be floats."

    if maxiter <= 0:
        return "Invalid input: maxiter must be positive."
    if popsize <= 0:
        return "Invalid input: popsize must be positive."
    if tol <= 0:
        return "Invalid input: tol must be positive."
    if not (0.0 < recombination <= 1.0):
        return "Invalid input: recombination must be in (0, 1]."

    # Handle mutation parameter: scalar or 2D list [[value]] or [[min, max]]
    mutation = to2d(mutation)
    if isinstance(mutation, list) and len(mutation) > 0 and isinstance(mutation[0], list):
        flat_mutation = mutation[0]
        if len(flat_mutation) == 1:
            try:
                mutation_value = float(flat_mutation[0])
            except (TypeError, ValueError):
                return "Invalid input: mutation must be numeric."
            if not (0 < mutation_value <= 2):
                return "Invalid input: mutation scalar must be in (0, 2]."
            mutation_param = mutation_value
        elif len(flat_mutation) == 2:
            try:
                mutation_low = float(flat_mutation[0])
                mutation_high = float(flat_mutation[1])
            except (TypeError, ValueError):
                return "Invalid input: mutation range must be numeric."
            if not (0 <= mutation_low <= mutation_high <= 2):
                return "Invalid input: mutation range must satisfy 0 <= low <= high <= 2."
            mutation_param = (mutation_low, mutation_high)
        else:
            return "Invalid input: mutation must be [[value]] or [[min, max]]."
    elif isinstance(mutation, (int, float)):
        try:
            mutation_value = float(mutation)
        except (TypeError, ValueError):
            return "Invalid input: mutation must be numeric."
        if not (0 < mutation_value <= 2):
            return "Invalid input: mutation scalar must be in (0, 2]."
        mutation_param = mutation_value
    else:
        return "Invalid input: mutation must be a float or 2D list [[value]] or [[min, max]]."

    rng_seed = None
    if seed is not None:
        try:
            rng_seed = int(seed)
        except (TypeError, ValueError):
            return "Invalid input: seed must be an integer."

    if not isinstance(polish, bool):
        return "Invalid input: polish must be a boolean."

    safe_globals = {
        "math": math,
        "np": np,
        "numpy": np,
        "__builtins__": {},
    }
    safe_globals.update({
        name: getattr(math, name)
        for name in dir(math)
        if not name.startswith("_")
    })
    safe_globals.update({
        "sin": np.sin,
        "cos": np.cos,
        "tan": np.tan,
        "asin": np.arcsin,
        "arcsin": np.arcsin,
        "acos": np.arccos,
        "arccos": np.arccos,
        "atan": np.arctan,
        "arctan": np.arctan,
        "sinh": np.sinh,
        "cosh": np.cosh,
        "tanh": np.tanh,
        "exp": np.exp,
        "log": np.log,
        "ln": np.log,
        "log10": np.log10,
        "sqrt": np.sqrt,
        "abs": np.abs,
        "pow": np.power,
        "pi": math.pi,
        "e": math.e,
    })

    x0_vector = [np.mean([b[0], b[1]]) for b in processed_bounds]

    try:
        initial_eval = eval(func_expr, safe_globals, {"x": x0_vector})
        float(initial_eval)
    except Exception as exc:
        return f"Error: Invalid expression at initial guess: {exc}"

    def objective(x_vector):
        try:
            local_x = [float(val) for val in np.atleast_1d(x_vector)]
            result = eval(func_expr, safe_globals, {"x": local_x})
            numeric_value = float(result)
        except Exception as exc:
            raise ValueError(f"Error evaluating func_expr: {exc}")
        if not math.isfinite(numeric_value):
            raise ValueError("Objective evaluation produced NaN or infinity.")
        return numeric_value

    try:
        result = scipy_differential_evolution(
            objective,
            bounds=processed_bounds,
            strategy=strategy,
            maxiter=maxiter,
            popsize=popsize,
            tol=tol,
            mutation=mutation_param,
            recombination=recombination,
            seed=rng_seed,
            polish=polish,
        )
    except ValueError as exc:
        return f"Error during diff_evolution: {exc}"
    except Exception as exc:
        return f"Error during diff_evolution: {exc}"

    if not result.success:
        return f"Error during diff_evolution: {result.message}"

    try:
        solution_vector = [float(val) for val in np.atleast_1d(result.x)]
    except (TypeError, ValueError):
        return "Error during diff_evolution: could not convert solution to floats."

    objective_value = float(result.fun)
    return [solution_vector + [objective_value]]

Overview

Excel Usage

Examples

Python Code

Online Calculator