如何在Python优化中考虑限制条件的情况下找到全局最小值？

Question

我有一个包含64个变量的Python函数，我尝试使用minimise函数中的L-BFGS-B方法进行优化，但是这种方法对初始猜测值有很强的依赖性，并且无法找到全局最小值。

但我喜欢它设置变量边界的能力。有没有一种方法/函数可以在具有变量边界的情况下找到全局最小值?

- dilyar

---

我怀疑http://math.stackexchange.com/是一个更适合询问这类问题的地方。 - 9000

您可以稍微描述一下您的函数 - 平滑/梯度/海森矩阵吗？如果您可以将其公式化为平方和，请参见 scipy-optimize-leastsq-with-bound-constraints 。另请参见 scicomp.stackexchange.com/search?q=bfgs。 - denis

我正在设计8条三维空间中的贝塞尔曲线，每条曲线有6个控制点。要最小化的函数是这些曲线的优点函数，它是从曲线中推导出的4个不同参数（长度、曲率半径、接近度、高度顺序）的线性组合。到目前为止，我已经尝试了scipy.minimize()和basinhopping，但仍然无法找到全局最小值。 - dilyar

Answer 1

这可以使用scipy.optimize.basinhopping来完成。Basinhopping是一个函数，旨在找到目标函数的全局最小值。它使用scipy.optimize.minimize函数进行重复最小化，并在每次最小化后在坐标空间中随机移动一步。通过使用实现边界的最小化器之一（例如L-BFGS-B），Basinhopping仍然可以尊重边界约束。以下是一些代码，演示如何实现此操作。

# an example function with multiple minima
def f(x): return x.dot(x) + sin(np.linalg.norm(x) * np.pi)

# the starting point
x0 = [10., 10.]

# the bounds
xmin = [1., 1.]
xmax = [11., 11.]

# rewrite the bounds in the way required by L-BFGS-B
bounds = [(low, high) for low, high in zip(xmin, xmax)]

# use method L-BFGS-B because the problem is smooth and bounded
minimizer_kwargs = dict(method="L-BFGS-B", bounds=bounds)
res = basinhopping(f, x0, minimizer_kwargs=minimizer_kwargs)
print res

上述代码可以用于简单的情况，但是如果 basinhopping 的随机位移例程将您带到了禁止区域，则仍然可能会出现问题。幸运的是，这可以通过使用关键字 take_step 传递自定义步骤取样例程来覆盖。

class RandomDisplacementBounds(object):
    """random displacement with bounds"""
    def __init__(self, xmin, xmax, stepsize=0.5):
        self.xmin = xmin
        self.xmax = xmax
        self.stepsize = stepsize

    def __call__(self, x):
        """take a random step but ensure the new position is within the bounds"""
        while True:
            # this could be done in a much more clever way, but it will work for example purposes
            xnew = x + np.random.uniform(-self.stepsize, self.stepsize, np.shape(x))
            if np.all(xnew < self.xmax) and np.all(xnew > self.xmin):
                break
        return xnew

# define the new step taking routine and pass it to basinhopping
take_step = RandomDisplacementBounds(xmin, xmax)
result = basinhopping(f, x0, niter=100, minimizer_kwargs=minimizer_kwargs,
                      take_step=take_step)
print result

- Jacob Stevenson

---

一种没有循环的变体：return np.clip( x + random(...), xmin, xmax )。（在所有维度中使用相同的步长，预缩放以使边界框大致为正方形。） - denis

你好。我认为在__call__中，你需要self.xmin和self.xmax而不是xmin和xmax？ - CPBL

Answer 2

针对任何优化器的调试和可视化，以下是一些常识性建议：

你的目标函数和约束是否合理？
如果目标函数是一个求和形式，比如 f() + g()， 请分别为所有的x在 "fx-opt.nptxt"（下面）中打印它们; 如果f()占总和的 99%，而 g() 仅占 1%，请进行调查。

约束条件：有多少个xfinal中的组成部分x_i被限制在边界上， 如x_i <= lo_i或>= hi_i？

---

你的函数在全局范围内有多么崎岖不平？
使用几个随机起点运行，并保存结果以进行分析/绘图：

title = "%s  n %d  ntermhess %d  nsample %d  seed %d" % (  # all params!
    __file__, n, ntermhess, nsample, seed )
print title
...
np.random.seed(seed)  # for reproducible runs
np.set_printoptions( threshold=100, edgeitems=10, linewidth=100,
        formatter = dict( float = lambda x: "%.3g" % x ))  # float arrays %.3g

lo, hi = bounds.T  # vecs of numbers or +- np.inf
print "lo:", lo
print "hi:", hi

fx = []  # accumulate all the final f, x
for jsample in range(nsample):
        # x0 uniformly random in box lo .. hi --
    x0 = lo + np.random.uniform( size=n ) * (hi - lo)

    x, f, d = fmin_l_bfgs_b( func, x0, approx_grad=1,
                m=ntermhess, factr=factr, pgtol=pgtol )
    print "f: %g  x: %s  x0: %s" % (f, x, x0)
    fx.append( np.r_[ f, x ])

fx = np.array(fx)  # nsample rows, 1 + dim cols
np.savetxt( "fx-opt.nptxt", fx, fmt="%8.3g", header=title )  # to analyze / plot

ffinal = fx[:,0]
xfinal = fx[:,1:]
print "final f values, sorted:", np.sort(ffinal)
jbest = ffinal.argmin()
print "best x:", xfinal[jbest]

如果一些ffinal值看起来相当不错， 尝试更多接近那些值的随机起点——这比纯随机要好。

如果x是曲线或其他实际对象，绘制最佳的几个x0和xfinal。
(一个经验法则是在d维度中使用nsample ~ 5*d或10*d。 过于缓慢或过多？减少maxiter/maxeval，减少ftol——对于这种探索，你不需要ftol1e-6。)

如果您想要可重复的结果， 那么您必须在title和派生文件和图表中列出所有相关参数。 否则，您会问“这是从哪里来的？”

---

您的函数在epsilon尺度为10^-6时有多崎岖？
有时，近似梯度的方法会返回它们的最后估计， 但如果没有：

from scipy.optimize._numdiff import approx_derivative  # 3-point, much better than
## from scipy.optimize import approx_fprime
for eps in [1e-3, 1e-6]:
    grad = approx_fprime( x, func, epsilon=eps )
    print "approx_fprime eps %g: %s" % (eps, grad)

如果优化器在停止之前梯度估计不够准确或者出现起伏，那么你是看不到的。那么你需要保存所有中间值[f, x, approx_fprime]来观察它们；在Python中很容易做到-如果这不清楚，请告诉我。

在某些问题领域，倒退并从一个假定的xmin重新开始是很常见的。例如，如果你在乡间小道上迷路了，首先找到一条主要道路，然后从那里重新开始。

总结： 不要指望任何黑匣子优化器在处理大规模波动或浅尺度波动或两者都有的函数时能够正常工作。 投资于测试脚手架和以可视方式观察优化器正在进行的操作。

- denis

Answer 3

自您提出此问题以来，全局优化方面已经有了一些不错的进展，这可能对您有所帮助。特别是，我想引起您对SHGO算法（package）的注意，它现在也是scipy.optimize的标准选项之一。然而，如果您真的无法减少搜索空间的维数，它可能会有些困难。

您可以尝试一些经典的模式搜索或PySOT中的代理方法，它们也表现良好。如果您真的卡住了，考虑使用optuna或者如果您非常绝望，可以使用hyperopt。

我目前的选择是DLIB和Nevergrad。两者都非常快速。如果您需要一个没有任何实际依赖关系的选项，可以看看freelunch。

- Peter Cotton

Answer 4

非常感谢您的详细回复，但由于我对Python还比较陌生，不太知道如何将代码应用到我的程序中。以下是我尝试进行优化的代码：

x0=np.array((10, 13, f*2.5, 0.08,    10, f*1.5,  0.06, 20, 
             10, 14, f*2.5, 0.08,    10, f*1.75, 0.07, 20,
             10, 15, f*2.5, 0.08,    10, f*2,    0.08, 20,
             10, 16, f*2.5, 0.08,    10, f*2.25, 0.09, 20,
             10, 17, f*2.5, -0.08,    10, f*2.5, -0.06, 20,
             10, 18, f*2.5, -0.08,    10, f*2.75,-0.07, 20,
             10, 19, f*2.5, -0.08,    10, f*3,   -0.08, 20,
             10, 20, f*2.5, -0.08,    10, f*3.25,-0.09, 20))

# boundary for each variable, each element in this restricts the corresponding element     above
bnds=((1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
  (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
  (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
  (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
  (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
  (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
  (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), 
  (1,12), (1,35), (0,f*6.75), (-0.1, 0.1),(1,35), (0,f*6.75), (-0.1, 0.1),(13, 35), )

from scipy.optimize import basinhopping
from scipy.optimize import minimize

merit=a*meritoflength + b*meritofROC + c*meritofproximity +d*(distancetoceiling+distancetofloor)+e*heightorder
minimizer_kwargs = {"method": "L-BFGS-B", "bounds": bnds, "tol":1e0}
ret = basinhopping(merit_function, x0, minimizer_kwargs=minimizer_kwargs, niter=10, T=0.01)

zoom = ret['x']
res = minimize(merit_function, zoom, method = 'L-BFGS-B', bounds=bnds, tol=1e-5)
print res

优点函数将x0与其他一些值结合起来形成8条曲线的6个控制点，然后计算它们的长度、曲率半径等。它将这些参数的线性组合与一些权重返回为最终的优点。

我使用低精度的basinhopping找到了一些极小值，然后使用minimize增加了最低极小值的精度。

附注：我正在运行的平台是Enthoght Canopy 1.3.0，numpy 1.8.0，scipy 0.13.2，mac 10.8.3。

- dilyar

---

你的函数“跳跃”穿过盆地的次数有多少？ 这对那些可能尝试此操作的人很有用。 还有一个一般性问题：控制点是否限制在一个大框中，还是每个点都在自己的1/6中？如果是一个大框，则解空间比需要的大720倍。 - denis

如何检查函数跳转次数？对于6个控制点，分别为P1、P2、P3、P4、P5、P6。在这些点中，P1、P6是固定的，P2y、P2z、P5y、p5z是固定的，因此变量是P2x、P3x、P3y、P3z、P4x、P4y、P4z、P5x。现在我相信P3和P4共享同一个大盒子，而P2、P5则沿着单一的x轴移动。 - dilyar

我建议提出一个新问题，比如“如何检查函数在scipy.optimize.basinhopping中的跳跃方式”。但是你已经得到了合理的结果，是吗？ - denis