scipy.optimize.

differential_evolution#

scipy.optimize.differential_evolution(func, bounds, args=(), strategy='best1bin', maxiter=1000, popsize=15, tol=0.01, mutation=(0.5, 1), recombination=0.7, rng=None, callback=None, disp=False, polish=True, init='latinhypercube', atol=0, updating='immediate', workers=1, constraints=(), x0=None, *, integrality=None, vectorized=False, seed=None)[source]#

寻找多元函数的全局最小值。

微分进化方法 [1] 本质上是随机的。它不使用梯度方法来寻找最小值,可以搜索候选空间的巨大区域,但通常比传统的基于梯度的技术需要更多的函数评估次数。

该算法由 Storn 和 Price [2] 提出。

参数:
funccallable

需要最小化的目标函数。必须采用 f(x, *args) 的形式,其中 x 是 1-D 数组形式的参数,args 是用于完全指定该函数所需的任何额外固定参数的元组。参数数量 N 等于 len(x)

boundssequence or Bounds

变量的边界。有两种方法指定边界

  1. Bounds 类的实例。

  2. (min, max) 对,分别对应 x 中的每个元素,定义了 func 优化参数的有限下界和上界。

边界的总数用于确定参数的数量 N。如果某些参数的边界相等,则自由参数的总数为 N - N_equal

args元组 (tuple),可选

完全指定目标函数所需的任何额外固定参数。

strategy{str, callable}, 可选

所使用的微分进化策略。应该是以下之一:

  • ‘best1bin’

  • ‘best1exp’

  • ‘rand1bin’

  • ‘rand1exp’

  • ‘rand2bin’

  • ‘rand2exp’

  • ‘randtobest1bin’

  • ‘randtobest1exp’

  • ‘currenttobest1bin’

  • ‘currenttobest1exp’

  • ‘best2exp’

  • ‘best2bin’

默认值为 ‘best1bin’。可实现的策略在“注释”中概述。或者,可以通过提供一个用于构建试验向量的可调用对象来定制微分进化策略。该可调用对象必须采用 strategy(candidate: int, population: np.ndarray, rng=None) 的形式,其中 candidate 是一个整数,指定正在演化的种群成员索引;population 是一个形状为 (S, N) 的数组,包含所有种群成员(其中 S 是种群总大小);rng 是求解器内部使用的随机数生成器。candidate 将在范围 [0, S) 内。strategy 必须返回一个形状为 (N,) 的试验向量。该试验向量的适应度将与 population[candidate] 的适应度进行比较。

版本 1.12.0 更改:支持通过可调用对象自定义进化策略。

maxiter整数,可选

整个种群进化的最大代数。函数评估的最大次数(不进行优化调整/polished)为:(maxiter + 1) * popsize * (N - N_equal)

popsizeint, 可选

用于设置总种群大小的乘数。种群拥有 popsize * (N - N_equal) 个个体。如果通过 init 关键字提供了初始种群,则此关键字将被覆盖。使用 init='sobol' 时,种群大小计算为 popsize * (N - N_equal) 之后的下一个 2 的幂。

tol浮点数,可选

收敛的相对容差。当 np.std(population_energies) <= atol + tol * np.abs(np.mean(population_energies)) 时求解停止,其中 atoltol 分别为绝对容差和相对容差。

mutationfloat 或 tuple(float, float), 可选

变异常数。在文献中也称为微分权重,记为 \(F\)。如果指定为浮点数,应在范围 [0, 2) 内。如果指定为元组 (min, max),则采用抖动(dithering)。抖动会随机改变每一代的变异常数。该代的变异常数取自 U[min, max)。抖动可以显著加速收敛。增加变异常数会扩大搜索半径,但会减慢收敛速度。

recombinationfloat, 可选

重组常数,应在范围 [0, 1] 内。在文献中也称为交叉概率,记为 CR。增加该值允许更多突变体进入下一代,但会增加种群不稳定的风险。

rng{None, int, numpy.random.Generator}, 可选

如果通过关键字传递 rng,则 numpy.random.Generator 以外的类型将被传递给 numpy.random.default_rng 以实例化 Generator。如果 rng 已经是 Generator 实例,则使用提供的实例。指定 rng 以获得可重复的函数行为。

如果此参数通过位置传递或 seed 通过关键字传递,则 seed 参数的旧行为适用

  • 如果 seed 为 None(或 numpy.random),则使用 numpy.random.RandomState 单例。

  • 如果 seed 是一个 int,则使用一个新的 RandomState 实例,并以 seed 进行播种。

  • 如果 seed 已经是 GeneratorRandomState 实例,则使用该实例。

版本 1.15.0 中更改: 作为从 numpy.random.RandomState 过渡到 numpy.random.GeneratorSPEC-007 计划的一部分,此关键字已从 seed 更改为 rng。在过渡期间,两个关键字都将继续有效,尽管每次只能指定其中一个。过渡期结束后,使用 seed 关键字的函数调用将发出警告。上述已概述了 seedrng 的行为,但在新代码中应仅使用 rng 关键字。

disp布尔值,可选

在每次迭代时打印评估的 func

callback可调用对象,可选

每次迭代后调用的可调用对象。具有以下签名:

callback(intermediate_result: OptimizeResult)

其中 intermediate_result 是一个关键字参数,包含一个具有 xfun 属性的 OptimizeResult,即目前找到的最佳解和目标函数值。请注意,参数名称必须是 intermediate_result,以便回调函数接收到 OptimizeResult

回调函数也支持如下签名:

callback(x, convergence: float=val)

val 表示种群收敛的分数值。当 val 大于 1.0 时,函数停止。

通过内省(Introspection)来确定调用哪种签名。

如果回调抛出 StopIteration 或返回 True,全局最小化将停止;任何优化调整(polishing)仍会执行。

版本 1.12.0 更改:回调接受 intermediate_result 关键字。

polish{bool, callable}, 可选

如果为 True(默认),则在最后使用带有 L-BFGS-B 方法的 scipy.optimize.minimize 来优化调整最佳种群成员,这可以轻微改善最小化效果。如果研究的是约束问题,则使用 trust-constr 方法。对于带有许多约束的大型问题,由于雅可比计算,优化调整可能需要很长时间。或者,提供一个具有类似 minimize 签名的可调用对象 polish_func(func, x0, **kwds) 并返回一个 OptimizeResult。这允许用户精细控制优化调整的发生方式。boundsconstraints 将出现在 kwds 中。可以使用 functools.partialpolish_func 提供额外关键字。用户有责任确保优化函数遵守边界、任何约束(包括整数约束),并且在 OptimizeResult 中设置了适当的属性,例如 funxnfevjac

版本 1.15.0 更改:如果指定了 workers,则会将封装 func 的映射类可调用对象提供给 minimize,而不是直接使用 func。这允许调用者控制调用实际运行的方式和位置。

版本 1.17.0 更改:可以提供一个遵守 minimize 签名的可调用对象来优化调整最佳种群成员。

initstr 或 array-like, 可选

指定执行哪种类型的种群初始化。应该是以下之一:

  • ‘latinhypercube’

  • ‘sobol’

  • ‘halton’

  • ‘random’

  • 指定初始种群的数组。数组形状应为 (S, N),其中 S 是种群总大小,N 是参数数量。

init 在使用前会截断至 bounds

默认值为 ‘latinhypercube’。拉丁超立方采样(Latin Hypercube sampling)试图最大化对可用参数空间的覆盖。

‘sobol’ 和 ‘halton’ 是更优的替代方案,能进一步最大化对参数空间的覆盖。‘sobol’ 将强制执行一个初始种群大小,计算为 popsize * (N - N_equal) 之后的下一个 2 的幂。‘halton’ 没有要求,但效率稍低。详情请参阅 scipy.stats.qmc

‘random’ 随机初始化种群 - 缺点是可能发生聚类,导致无法覆盖整个参数空间。使用数组指定种群可以用于例如在已知解存在的位置创建一个紧密的初始猜测簇,从而减少收敛时间。

atolfloat, 可选

收敛的绝对容差,当 np.std(population_energies) <= atol + tol * np.abs(np.mean(population_energies)) 时求解停止,其中 atoltol 分别为绝对容差和相对容差。

updating{‘immediate’, ‘deferred’}, 可选

如果为 'immediate',则在单代内持续更新最佳解向量 [4]。这可能导致更快的收敛,因为试验向量可以利用最佳解的持续改进。使用 'deferred',最佳解向量每代更新一次。只有 'deferred' 与并行化或向量化兼容,workersvectorized 关键字可以覆盖此选项。

在 1.2.0 版本中添加。

workersint or map-like callable, optional

如果 workers 是 int,则种群被划分为 workers 个部分并行评估(使用 multiprocessing.Pool)。提供 -1 可使用所有可用的 CPU 核心。或者,提供一个映射类(map-like)可调用对象,例如 multiprocessing.Pool.map,用于并行评估种群。评估执行方式为 workers(func, iterable)。如果 workers != 1,此选项会将 updating 关键字覆盖为 updating='deferred'。如果 workers != 1,此选项会覆盖 vectorized 关键字。要求 func 可被 pickle 处理。

在 1.2.0 版本中添加。

constraints{NonLinearConstraint, LinearConstraint, Bounds}

求解器的约束,在 bounds 关键字应用的基础之上。使用 Lampinen 的方法 [5]

版本 1.4.0 中新增。

x0None 或 array-like, 可选

提供最小化的初始猜测。一旦种群初始化,此向量将替换第一个(最佳)成员。即使通过 init 给出了初始种群,也会执行此替换。x0.shape == (N,)

版本 1.7.0 中新增。

integrality1-D array, 可选

对于每个决策变量,一个布尔值,指示该决策变量是否被约束为整数值。该数组会广播至 (N,)。如果任何决策变量被约束为整数,它们在优化调整期间将不会更改。仅使用位于下界和上界之间的整数值。如果界限之间没有整数值,则会引发 ValueError

版本 1.9.0 中新增。

vectorized布尔值,可选

如果 vectorized is Truefunc 将接收一个形状为 x.shape == (N, S)x 数组,并预期返回一个形状为 (S,) 的数组,其中 S 是要计算的解向量数量。如果应用了约束,每个用于构建 Constraint 对象的函数都应接受形状为 x.shape == (N, S)x 数组,并返回一个形状为 (M, S) 的数组,其中 M 是约束分量的数量。此选项是 workers 提供的并行化方案的一种替代方案,可以通过减少多次函数调用带来的解释器开销来提高优化速度。如果 workers != 1,则忽略此关键字。此选项会将 updating 关键字覆盖为 updating='deferred'。请参阅注释部分,进一步讨论何时使用 'vectorized' 以及何时使用 'workers'

版本 1.9.0 中新增。

返回:
resOptimizeResult

OptimizeResult 对象表示的优化结果。重要的属性有:x(解数组)、success(布尔标志,指示优化器是否成功退出)、message(描述终止原因)、population(存在于种群中的解向量)以及 population_energiespopulation 中每个条目的目标函数值)。有关其他属性的描述,请参阅 OptimizeResult。如果采用了 polish 且通过优化调整获得了更低的最小值,则 OptimizeResult 还会包含 jac 属性。如果最终解不满足应用的约束,success 将为 False

附注

微分进化是一种基于种群的随机方法,适用于全局优化问题。在每次通过种群时,该算法通过与其他候选解混合来变异每个候选解,从而创建一个试验候选者。有几种用于创建试验候选者的策略 [3],某些策略比其他策略更适合某些问题。‘best1bin’ 策略是许多系统的良好起点。在此策略中,随机选择种群中的两个成员。它们的差值用于变异目前最好的成员(即 ‘best1bin’ 中的 ‘best’)\(x_0\)

\[b' = x_0 + F \cdot (x_{r_0} - x_{r_1})\]

其中 \(F\)mutation 参数。然后构建一个试验向量。从随机选择的第 i 个参数开始,试验向量被依次(以模数方式)填充来自 b' 或原始候选者的参数。是否使用 b' 或原始候选者的选择是通过二项分布(‘best1bin’ 中的 ‘bin’)做出的——生成一个 [0, 1) 范围内的随机数。如果该数小于 recombination 常数,则从 b' 加载参数,否则从原始候选者加载。始终从 b' 加载一个随机选择的参数。对于二项交叉,这是一个单一的随机参数。对于指数交叉,这是来自 b' 的连续参数序列的起点。一旦构建好试验候选者,就会评估其适应度。如果试验候选者优于原始候选者,它就会取代原来的位置。如果它也优于整体最佳候选者,它也会取代那个位置。

其他可用策略在 Qiang 和 Mitchell (2014) [3] 中概述。

  • rand1 : \(b' = x_{r_0} + F \cdot (x_{r_1} - x_{r_2})\)

  • rand2 : \(b' = x_{r_0} + F \cdot (x_{r_1} + x_{r_2} - x_{r_3} - x_{r_4})\)

  • best1 : \(b' = x_0 + F \cdot (x_{r_0} - x_{r_1})\)

  • best2 : \(b' = x_0 + F \cdot (x_{r_0} + x_{r_1} - x_{r_2} - x_{r_3})\)

  • currenttobest1 : \(b' = x_i + F \cdot (x_0 - x_i + x_{r_0} - x_{r_1})\)

  • randtobest1 : \(b' = x_{r_0} + F \cdot (x_0 - x_{r_0} + x_{r_1} - x_{r_2})\)

其中整数 \(r_0, r_1, r_2, r_3, r_4\) 从区间 [0, NP) 中随机选择,NP 为总种群大小,原始候选者索引为 i。用户可以通过向 strategy 提供可调用对象来完全自定义试验候选者的生成。

为了提高找到全局最小值的机会,请使用更高的 popsize 值,结合更高的 mutation 和(抖动),以及更低的 recombination 值。这会起到扩大搜索半径但减慢收敛速度的效果。

默认情况下,最佳解向量在单次迭代中持续更新(updating='immediate')。这是原始微分进化算法的一种修改 [4],由于试验向量可以立即受益于改进的解,这可能导致更快的收敛。要使用原始的 Storn 和 Price 行为(每迭代更新一次最佳解),请设置 updating='deferred''deferred' 方法与并行化和向量化(workersvectorized 关键字)都兼容。这些方法可以通过更有效地使用计算机资源来提高最小化速度。'workers' 将计算分发到多个处理器。默认使用 Python multiprocessing 模块,但也可能使用其他方法,例如在集群上使用的消息传递接口 (MPI) [6] [7]。这些方法的开销(创建新进程等)可能是巨大的,这意味着计算速度不一定会随着所用处理器数量的增加而线性扩展。并行化最适合计算昂贵的目标函数。如果目标函数成本较低,则 'vectorized' 可能会通过在每次迭代中仅调用一次目标函数(而不是为所有种群成员多次调用)来提供帮助;减少了解释器开销。

版本 0.15.0 中新增。

参考文献

[1]

Differential evolution, Wikipedia, http://en.wikipedia.org/wiki/Differential_evolution

[2]

Storn, R and Price, K, Differential Evolution - a Simple and Efficient Heuristic for Global Optimization over Continuous Spaces, Journal of Global Optimization, 1997, 11, 341 - 359.

[3] (1,2)

Qiang, J., Mitchell, C., A Unified Differential Evolution Algorithm for Global Optimization, 2014, https://www.osti.gov/servlets/purl/1163659

[4] (1,2)

Wormington, M., Panaccione, C., Matney, K. M., Bowen, D. K., - Characterization of structures from X-ray scattering data using genetic algorithms, Phil. Trans. R. Soc. Lond. A, 1999, 357, 2827-2848

[5]

Lampinen, J., A constraint handling approach for the differential evolution algorithm. Proceedings of the 2002 Congress on Evolutionary Computation. CEC’02 (Cat. No. 02TH8600). Vol. 2. IEEE, 2002.

示例

让我们考虑最小化 Rosenbrock 函数的问题。此函数在 scipy.optimize 中的 rosen 中实现。

>>> import numpy as np
>>> from scipy.optimize import rosen, differential_evolution
>>> bounds = [(0,2), (0, 2), (0, 2), (0, 2), (0, 2)]
>>> result = differential_evolution(rosen, bounds)
>>> result.x, result.fun
(array([1., 1., 1., 1., 1.]), 1.9216496320061384e-19)

现在重复上述过程,但使用并行化。

>>> result = differential_evolution(rosen, bounds, updating='deferred',
...                                 workers=2)
>>> result.x, result.fun
(array([1., 1., 1., 1., 1.]), 1.9216496320061384e-19)

让我们执行一个受约束的最小化。

>>> from scipy.optimize import LinearConstraint, Bounds

我们添加一个约束,即 x[0]x[1] 的和必须小于或等于 1.9。这是一个线性约束,可以写为 A @ x <= 1.9,其中 A = array([[1, 1]])。这可以编码为 LinearConstraint 实例。

>>> lc = LinearConstraint([[1, 1]], -np.inf, 1.9)

使用 Bounds 对象指定界限。

>>> bounds = Bounds([0., 0.], [2., 2.])
>>> result = differential_evolution(rosen, bounds, constraints=lc,
...                                 rng=1)
>>> result.x, result.fun
(array([0.96632622, 0.93367155]), 0.0011352416852625719)

接下来寻找 Ackley 函数的最小值(https://en.wikipedia.org/wiki/Test_functions_for_optimization)。

>>> def ackley(x):
...     arg1 = -0.2 * np.sqrt(0.5 * (x[0] ** 2 + x[1] ** 2))
...     arg2 = 0.5 * (np.cos(2. * np.pi * x[0]) + np.cos(2. * np.pi * x[1]))
...     return -20. * np.exp(arg1) - np.exp(arg2) + 20. + np.e
>>> bounds = [(-5, 5), (-5, 5)]
>>> result = differential_evolution(ackley, bounds, rng=1)
>>> result.x, result.fun
(array([0., 0.]), 4.440892098500626e-16)

Ackley 函数以向量化方式编写,因此可以使用 'vectorized' 关键字。请注意函数评估次数的减少。

>>> result = differential_evolution(
...     ackley, bounds, vectorized=True, updating='deferred', rng=1
... )
>>> result.x, result.fun
(array([0., 0.]), 4.440892098500626e-16)

最终的优化调整步骤可以通过提供模仿 minimize 接口的可调用对象来定制。用户有责任确保最小化器遵守任何边界和约束。

>>> from functools import partial
>>> from scipy.optimize import minimize
>>> # supply extra parameters to the polishing function using partial
>>> polish_func = partial(minimize, method="SLSQP")
>>> result = differential_evolution(
...     ackley, bounds, vectorized=True, updating='deferred', rng=1,
...     polish=polish_func
... )
>>> result.x, result.fun
(array([0., 0.]), 4.440892098500626e-16)

以下自定义策略函数模仿 ‘best1bin’。

>>> def custom_strategy_fn(candidate, population, rng=None):
...     parameter_count = population.shape[-1]
...     mutation, recombination = 0.7, 0.9
...     trial = np.copy(population[candidate])
...     fill_point = rng.choice(parameter_count)
...
...     pool = np.arange(len(population))
...     rng.shuffle(pool)
...
...     # two unique random numbers that aren't the same, and
...     # aren't equal to candidate.
...     idxs = []
...     while len(idxs) < 2 and len(pool) > 0:
...         idx = pool[0]
...         pool = pool[1:]
...         if idx != candidate:
...             idxs.append(idx)
...
...     r0, r1 = idxs[:2]
...
...     bprime = (population[0] + mutation *
...               (population[r0] - population[r1]))
...
...     crossovers = rng.uniform(size=parameter_count)
...     crossovers = crossovers < recombination
...     crossovers[fill_point] = True
...     trial = np.where(crossovers, bprime, trial)
...     return trial