使用 scipy.optimize.minimize 编写目标函数故障排除

Question

我正在尝试在参数化设计的背景下解决应用于桥梁的优化问题。目标函数是尝试将预制障碍物分割成最少量的独特障碍物。护栏有普通护栏和灯杆护栏两种。我分享了一张 3 跨桥的图片。

参数如下：

灯杆屏障宽度等于w

灯杆屏障按参数 s

均匀分布

灯杆起点由变量a定义。

桥的跨度，[l1,l2,...ln]

x 是所有方程都必须被整除的宽度。

这个想法是使用尽可能多的宽度为 x 的障碍和一些剩余长度来填充保持最小的跨度。我必须获得方程组 n 中的所有距离 [a,b,c,...f]（参见图片），其中我还在特殊障碍物 (s-w) 之间附加了光，并且有它们都可以被x整除。我还在 elements arr 数组中添加了常量，以防剩余长度不能被 x 整除，因此我们将用不同的屏障进行补偿。 （基本上，我们使用该常数来减去

这是我正在谈论的内容的草图： .

from scipy.optimize import minimize
import numpy as np

s = 12500  
spans = np.array([25450, 25100, 25450])
ejs = np.zeros(len(spans)-1)  # expansion or separation joints . ignore for now

min_bw, max_bw  = 0, 1200  #light pole barrier width

consts = abs((len(spans) * 2 - 1))  # Number of distance elements [b,c,d...]
init = np.array([1500, 1500, 6250]) #x, w, a + w/2

n = (spans - s / 2) // s  # consecutive spacings of barriers for a given span. At times we have more 
m = np.ones(len(spans) * 2)  # the length of the equation array that has to be divisible

def objective(arr, n, s, spans, ejs):

    x = arr[0]
    w = arr[1]
    a =  arr[2] # (the offset is until the mid axis of the light pole)
    m = np.ones(len(spans) * 2 + 1) # We make room for the last constraint h=(s-w)
    m[0] = a - w / 2 

    for i in range(1, len(spans) * 2):
        j = (i - 1) // 2
        if i % 2 == 0:
            m[i] = (s - w) - (m[i - 1] + ejs[j] + arr[2 + i])
        else:
            m[i] = spans[j] - (m[i - 1] + n[j] * s + w + arr[2 + i])
    m[-1] = s - w  # We also need the light distance between two poles to be divisible by x
    #print(m)
    return np.sum(m%x)

#Bounds
bounds = np.array([(1500, 2100), (1500, 2600), (3000, s)])
extras = np.repeat([(min_bw, max_bw)], consts, axis=0) # Add the bounds for each additional element
all_bounds = np.concatenate((bounds, extras), axis=0)


initial = np.pad(init, (0, consts), 'constant', constant_values=min_bw)

sol = minimize(fun=objective, x0=initial, bounds=all_bounds, method='Nelder-Mead', args=(n, s, spans, ejs))
print(np.round(sol.x))

#m = [ 5500. , 5950.,  5050.  ,6050.  ,4950. , 6500. ,11000.] #the print m statment

我似乎没有弄清楚的问题是 np.sum(m%x) 不返回 0！与绘图相比，我在 m 数组中得到了正确的距离。未满足整除方程约束，但函数表示已成功终止。我将不胜感激任何见解！非常感谢！

Answer 1

您遇到了两个问题。

1. m % x

   有一个令优化器感到困惑的导数。
2. 您正在解决的问题有很多局部最小值。

让我们首先为优化器提供有关如何优化函数的更多信息。

令人困惑的导数

函数

minimize()

是一个局部优化器。如果还有另一个最低限度，它必须爬上一座小山才能到达，它就无法到达。

假设您想要最小化函数 f(x) = x % 100，其中 x 介于 50 和 150 之间。这是该函数的图表。

假设优化器从 x = 90 开始。它查看导数，然后沿着斜率向左移动。然后，它陷入 x = 50。同样的事情也发生在你的问题中。

解决这个问题你可以做的是创建一个函数，奖励优化器接近 x = 100。我想出了这个函数，它是一个余弦波，当

variable

接近零时，它会接近零。

target

的倍数。

def modular_difference(variable, target):
    # Move valley very slightly forward
    overshoot_distances_amount = 1e-8
    distance_to_next = variable % target / target
    distance_to_next -= overshoot_distances_amount
    ret = -np.cos(2*np.pi*distance_to_next) + 1
    return ret

（注：

overshoot_distances_amount

的用途稍后解释。）

这给出了下图：

这告诉优化器，如果 x=99，它应该向 x=100 移动。

我通过以下方式将这个新功能集成到您的代码中。

首先，我更改了原始目标函数，使其同时返回 m 和 x。称之为

objective_inner()

。

def objective_inner(arr, n, s, spans, ejs):
    x = arr[0]
    w = arr[1]
    a =  arr[2] # (the offset is until the mid axis of the light pole)
    m = np.ones(len(spans) * 2 + 1) # We make room for the last constraint h=(s-w)
    m[0] = a - w / 2 

    for i in range(1, len(spans) * 2):
        j = (i - 1) // 2
        if i % 2 == 0:
            m[i] = (s - w) - (m[i - 1] + ejs[j] + arr[2 + i])
        else:
            m[i] = spans[j] - (m[i - 1] + n[j] * s + w + arr[2 + i])
    m[-1] = s - w  # We also need the light distance between two poles to be divisible by x
    return m, x

然后，我写了两个版本的目标函数。

def objective_differentiable(arr, n, s, spans, ejs):
    m, x = objective_inner(arr, n, s, spans, ejs)
    return modular_difference(m, x).sum()


def objective_old(arr, n, s, spans, ejs):
    m, x = objective_inner(arr, n, s, spans, ejs)
    return np.sum(m%x)

第一个函数是目标的一个版本，它具有很好的导数，使用正弦波。第二个功能与您最初的目标相同。

但是，仅修复导数还不够。它仍然陷入局部最小值。

全局优化器

某些函数无法通过局部优化器进行优化，因为优化器会陷入局部最小值。这是其中之一。

为了解决这个问题，我使用了 basinhopping 函数，该函数使用本地优化器，但也在随机方向上采取大步，希望找到更好的最小值。

我发现以下方法效果很好。

minimizer_kwargs = dict(bounds=all_bounds, method='Nelder-Mead', args=(n, s, spans, ejs))
sol = basinhopping(objective_differentiable, x0=initial, minimizer_kwargs=minimizer_kwargs, niter=100)
print("basinhopping sol", sol)
sol = minimize(fun=objective_old, x0=sol.x, bounds=all_bounds, method='Nelder-Mead', args=(n, s, spans, ejs))
print("minimize true sol", sol)

这会执行以下操作：

• 首先，它使用 Basinhopping 和 Nelder-Mead 作为本地优化器，优化目标的可微分版本。
• 然后，它将找到的解决方案传递给您最初的目标。我这样做是为了很容易比较它是否找到了比原来更好的解决方案。 （注意：这就是上面代码中存在

  overshoot_distances_amount

  的原因。这会选择一个稍微位于模数环绕点之后的解决方案，这使得优化器更容易找到

  objective_old()

  的精确最小值。）

这可以找到在您的原始目标上得分为 1e-4 或更低的解决方案。

这是它找到的解决方案之一：

array([1501.43791458, 1989.93457535, 5499.28104559,  449.99999353,
          0.000001  ,  100.00001273,    0.00000361,  450.00000564])