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.
Example 1: Sphere function with explicit 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 range
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 | 2.46519e-31 |
Example 3: Absolute value function without local 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 shifted quadratic 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 | 2.46519e-31 |
Python Code
Show 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, rand2bin, rand2exp, randtobest1bin, randtobest1exp, currenttobest1bin, 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.
"""
try:
def to2d(x):
return [[x]] if not isinstance(x, list) else x
if not isinstance(func_expr, str):
return "Error: func_expr must be a string."
if func_expr.strip() == "":
return "Error: 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 "Error: func_expr must reference variable x (e.g., x[0])."
bounds = to2d(bounds)
if not isinstance(bounds, list) or len(bounds) == 0:
return "Error: 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 "Error: each bound must be a [min, max] pair."
try:
lower = float(bound[0])
upper = float(bound[1])
except (TypeError, ValueError):
return "Error: bounds must contain numeric values."
if lower > upper:
return f"Error: lower bound must not exceed upper bound for variable index {idx}."
processed_bounds.append((lower, upper))
valid_strategies = {
'best1bin', 'best1exp', 'rand1exp', 'randtobest1exp', 'currenttobest1exp',
'best2exp', 'rand2exp', 'randtobest1bin', 'best2bin', 'rand2bin',
'rand1bin', 'currenttobest1bin'
}
if strategy not in valid_strategies:
return f"Error: 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 "Error: maxiter, popsize must be integers; tol, recombination must be floats."
if maxiter <= 0:
return "Error: maxiter must be positive."
if popsize <= 0:
return "Error: popsize must be positive."
if tol <= 0:
return "Error: tol must be positive."
if not (0.0 < recombination <= 1.0):
return "Error: recombination must be in (0, 1]."
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 "Error: mutation must be numeric."
if not (0 < mutation_value <= 2):
return "Error: 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 "Error: mutation range must be numeric."
if not (0 <= mutation_low <= mutation_high <= 2):
return "Error: mutation range must satisfy 0 <= low <= high <= 2."
mutation_param = (mutation_low, mutation_high)
else:
return "Error: mutation must be [[value]] or [[min, max]]."
elif isinstance(mutation, (int, float)):
try:
mutation_value = float(mutation)
except (TypeError, ValueError):
return "Error: mutation must be numeric."
if not (0 < mutation_value <= 2):
return "Error: mutation scalar must be in (0, 2]."
mutation_param = mutation_value
else:
return "Error: 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 "Error: seed must be an integer."
if not isinstance(polish, bool):
return "Error: 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]]
except Exception as e:
return f"Error: {str(e)}"